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ABSTRACT 


Seven algorithms for failure detection, isolation, and correction 
of redundant inertial instruments in the strqsdown dodecahedron 
configuration are competitively evaluated in a digital computer 
simulation that subjects them to identical environments. Their 
performance is compared in terms of orientation and inertial veloc- 
ity errors and in terms of missed and false alarms. 

The algorithms appear in the simulation program in modular form, 
so that they may be readily extracted for use elsewhere. The 
simulation program and its inputs and outputs are described. 

The algorithms, along with an eighth algorithm that was notsimu- 
lated, are also compared analytically to show the relationships 
among them. 


111 



PRECEDING PAGE BLANK NOT FILMED' 


TABLE OF CONTENTS 

Page 

1. INTRODUCTION 1 

1. 1 Background 1 

1. 2 What is a Failure? 2 

1. 3 Reliability and Performance 3 

1. 3. 1 Reliability with Perfect FDIC 5 

1. 3. 2 Performance with Perfect FDIC 6 

1. 4 The Algorithms 8 

1 . 5 The Simulation 1 0 

2. CONCLUSIONS 14 

3. CANDIDATE FDIC ALGORITHMS 16 

3. 1 Analytical Comparison of Algorithms 16 

3. 1. 1 General Considerations 16 

3. 1. 2 Derivation of Equations 20 

3. 1. 2. 1 Direct Test Signals 23 

3. 1.2.2 Variance Estimators 25 

3. 1. 2. 3 Indirect, Test Signals 29 

3. 1.2.4 Sequential Algorithm 32 

3. 1. 3 Multiple Simultaneous Failures 37 

3. 1. 4 Linearity Versus Nonlinearity 38 

3. 1. 5 V Space Logic 39 

3. 1. 6 Filtering 40 

3. 2 Detailed Algorithm Definitions 50 

3. 2. 1 Adaptive 66 Algorithm 50 

3. 2. 2 Fifteen Threshold Algorithm 52 

3. 2. 3 Squared Error Algorithm 57 

3. 2. 4 Bayesian Decision Theory Algorithm 63 

3. 2. 5 Maximum Likelihood Algorithm 69 

3. 2. 6 Minimax Algorithm 71 

3. 2. 7 Adaptive 72 Algorithm 77 

3. 2. 8 Sequential Algorithm 7 8 


Preceding page blank 


v 



TABLE OF CONTENTS (Continued) 


5. 


4. 1 Instrument Configuration 

4. 2 Unfailed Instrument Errors 

4. 2. 1 Gyro Errors 

4. 2. 2 Accelerometer Errors 

4. 3 Strapdown Algorithm 
4.4 Trajectory 

4.4. 1 Nominal Trajectory 
4. 4. 2 Vibrational Motions 
4. 5 Use of Inputs from Real Instruments 
FAILURE MODES 


5. 1 Gyro Failures 

5. 1. 1 Additive Failure Modes 

5. 1.2 Substitutional Failure Modes 

5. 2 Accelerometer Failures 

5. 2. 1 Additive Failure Modes 
5.2.2 Substitutional Failure Modes 


FAILSIM 

6. 1 Program Structure 

6. 2 Inputs 

6. 2. 1 Program Control 

6. 2. 2 Vehicle Vibrational Motions 

6. 2. 3 Instrument Orientation 
6. 2. 4 Unfailed Errors 
6. 2. 5 Failures 

6. 2. 6 External Strapdown Package 
6. 2. 7 Algorithms 
6.2.8 Nominal Trajectory 
6. 3 Outputs 

6. 3. 1 Inputs 

6. 3. 2 Events 

6. 3. 3 Figures of Merit 


Page 

r» o 
O O 

83 

86 

86 

87 

88 
89 
89 
89 
95 
98 
98 
98 

100 

100 

100 

101 

102 

102 

106 

106 

107 

107 

108 
109 
112 
113 
118 
1 19 
119 
1 19 
121 


vi 



TABLE OF CONTENTS (Continued) 


Page 

6. 3.4 Algorithm Execution Times 122 

6. 3. 5 Statistics 122 

6. 3. 6 Plot Tape 122 

7. SETTING THE PARAMETERS 124 

8. SIMULATION RESULTS 127 

8. 1 100-Sec Time Constant Runs 127 

8. 1. 1 First Failures 128 

8. 1.2 Second Failures 128 

8.2 10-Sec Time Constant Runs 128 

8.2. 1 First Failures 128 

8.2.2 Second Failures 128 

8.3 Double Simultaneous Failure Runs 132 

8. 4 Minor Cycle Time 134 

9. ALGORITHM SIZING AND TIMING 136 

9. 1 Algorithm Sizing 136 

9. 2 Algorithm Timing 136 

10. AREAS FOR FUTURE STUDY 138 

REFERENCES 140 

APPENDIX A. TYPICAL RUN 143 

APPENDIX B. THE SEQUENTIAL ALGORITHM IN STEADY STATE 156 


Vll 



ILLUSTRATIONS 


Page 

i_i FDIC Simulation Program (FAILSIM) 

Overall System Block Diagram 11 

3-1 Dodecahedron Instrument Configuration 22 

3-2 Low- Pass Prefilter Error Versus 

Failure Magnitude 42 

3-3 Gyro or Accelerometer Model 42 

3-4 Normalized Drift Rate Versus Normalized 

Time Constant 44 

3-5 Squared Error Algorithm Prefilter 45 

3-6 Prefilter Responses to a Step of Magnitude a 46 

3-7 Squared Error Algorithm Prefilter Error 

Versus Failure Magnitude 47 

3-8 Detection System for Positive Polarity 48 

3-9 Adaptive 72 Algorithm Filter Error 

Versus Failure Magnitude 49 

3-10 Adaptive 66 Algorithm Flow Diagram 51 

3-11 Fifteen Threshold Algorithm Flow Diagram 53 

3-12 Squared Error Algorithm Flow Diagram 61 

3-13 Bayesian Decision Theory Algorithm Flow Diagram 67 

3-14 Maximum Likelihood Algorithm Flow Diagram 70 

3-15 Minimax Algorithm Flow Diagram 73 

3- 16 Sequential Algorithm Flow Diagram 79 

4- 1 Vibrational Motion Generator 90 

4-2 Sampled Data Vibrational Motion Generator 94 

6-1 FAILSIM Flow Diagram 103 

6-2 Event Logic 121 

8-1 Orientation Error Magnitude Versus Failure 

Magnitude - First Gyro Failure, 

100-Sec Time Constant 129 

8-2 Orientation Error Magnitude Versus Failure . 

Magnitude - Second Gyro Failure, 

100-Sec Time Constant 129 

8-3 Orientation Error Magnitude Versus Failure 

Magnitude - First Gyro Failure, 

10-Sec Time Constant 130 

8-4 Orientation Error Magnitude Versus Failure 

Magnitude - Second Gyro Failure, 

10-Sec Time Constant 130 

viii 



ILLUSTRATIONS (Continued) 


8-5 Theoretical Error Magnitude Versus 

Failure Magnitude 


Page 

132 


IX 



TABLES 


Pag e 


1-1 Errors With Perfect FDIC 7 

4-1 Simulation - SIRU Correspondence 85 

7- 1 Adjustable Parameters in the Several Algorithms 124 

8- 1 FDIC Parameters 127 

8 -II FDIC Parameters 131 

8 -III FDIC Parameters 133 

8-IV Results of Double Simultaneous Failures 133 

8- V Summary of Events 134 

9- 1 Algorithm Sizing 136 

9-II Total Processor Time 137 


x 



ACKNOW LEDGEMENT 


This study would not have been possible without the 
assistance of the many people who provided me with 
information about the various algorithms, including 
T.-T. Chien, J. C. Deckert, R. E. Eckelkamp, J. T, 
Ephgrave, J. P. Gilmore, S. W. Gully, M. E. Jones, 
R. A, McKern, andR. T. Savely. However, I bear 
sole responsibility for the form the various algorithms 
take herein, especially since I found it necessary to 
modify them to some extent in order to make the 
competition fair. 


J. C. W. 


xi 



1. INTRODUCTION 


1. 1 BACKGROUND 

The notion of improving the reliability of a system by using 
redundant elements and some method of reorganization after a failure 
has been around for a long time. Since inertial instruments have never 
been as reliable as one might like, inertial navigators and attitude 
references have been prime candidates for this treatment. System level 
redundancy has been most common, with duplicated or triplicated gim- 
balled platforms or strapdown packages. Duplication permits automatic 
failure detection but requires some kind of external information for 
failure isolation. Triplication permits automatic failure detection and 
isolation by a simple majority voting scheme. However, a second fail- 
ure cannot always be isolated with triplicated systems. 

1 2 

Weiss and Nathan ' seem to have been the first to point out 
that, if six inertial instruments are arranged so that no three of their 
input axes are coplanar*, then first and second failures can be detected 
and isolated. Thus greater reliability can be achieved than with nine 
instruments arranged three per orthogonal axis. A third failure cannot 
be isolated without external information. 

3 4 5 6 

Ephgrave and Gilmore * ' have shown that the optimal 

arrangement for the six instruments is with their input axes perpendi- 
cular to the faces of a regular dodecahedron. Both gyros and accelero- 
meters can be arranged in this manner. The symmetry of the dodecahe- 
dron configuration maximizes accuracy for the worst cases of operation 
with a subset of the instruments. 

The simple majority voting scheme is no longer usable because 
the instrument outputs are not directly comparable, since their input 


It appears to be most practicable to embody this concept in a strapdown 
package, thus avoiding the problem of providing gimbal redundancy. 
However, it is possible, for example, to put three of the instruments on 
each of two gimballed platforms and use data crossfeeding and gimbal 
slaving techniques to enhance reliability. Only the strapdown case is 
considered here. 
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axes all point in different directions. Therefore, an algorithm must be 
devised to perform failure detection, isolation, and correction (FDIC), 
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extent to which the different algorithms have been reduced to practice 
varies widely from one to another. The purpose of this study is to com- 
pare an appropriate selection of FDIC algorithms with one another ana- 
lytically and by simulation to determine which is the best and to obtain 
insight into their modes of operation with the goal of combining their best 
features and remedying their shortcomings. 

The problem has two aspects. The first is the problem of detecting 
signals (the errors of the failed instruments) in the presence of noise 
(the errors of the unfailed instruments) and making a decision as to the 
existence of one or more failures. The second aspect is the problem of 
isolating the failed instruments. This problem is not trivial, because 
the available information is not always unambiguous in the case of two 
failures. Of course, simultaneous failures should be quite unlikely in 
a well-designed system, but one should not overlook the possibility that 
a second failure could occur during the finite time required to detect and 
isolate a first failure. It is also conceivable that a first failure might 
be small enough so as to go undetected and yet be able to interfere with 
the detection of a subsequent failure because of the ambiguity mentioned 
above. 

1. 2 WHAT IS A FAILURE? 

A definition of a failure might be: 

"A failure is an event wherein one or more com- 
ponents of an inertial instrument or its associated 
electronics ceases to function properly." 

This may be what intuition tells us a failure is, but it is not very useful 
in practice because we cannot, in general, tell whether or not each of the 
components of an instrument is functioning properly at every instant of 
time. 

A more operational definition might be: 

"A failure is an event -wherein the error in the 
output of an inertial instrument exceeds some 
predetermined specification. " 
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Unfortunately, the normal operation of an unfailed instrument may 
occasionally cause the specification to be exceeded during staging tran- 
sients, for example. Also it is mathematically convenient to model the 
errors of unfailed instruments with gaussian probability densities, which 
means that any specification has a finite probability of being exceeded 
without a failure. 

We shall cut the Gordian knot by defining a failure as: 

"A failure is an event wherein the failure simu- 
lator in the simulation program changes the state 
of an instrument from unfailed to failed with an 
accompanying change in the instrument error. " 

This definition makes up in convenience for what it lacks in 
elegance. 

1. 3 RELIABILITY AND PERFORMANCE 

In general, we consider reliability to be the probability that the 
system completes its mission without system failure (in some sense) 
and performance to be some measure of the errors of the unfailed sys* 
terns. The errors of the failed systems are large and ought not to be 
counted against performance. 

The nature of the interaction between reliability and performance 
depends upon the definition of reliability. A definition of reliability might 
be: 

"Reliability is the probability that less than three 
gyro failures and less than three accelerometer 
failures occur during the mission. " 

With this definition, the reliability is independent of the choice of algo- 
rithm. The relative merits of the different algorithms are then revealed 
by their performance, both in terms of missed and false alarms and in 
terms of the system velocity and orientation errors. 

Another definition of reliability might be: 

"Reliability is the probability that the system 
errors do not exceed their specifications dur- 
ing the mission. " 
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With this definition, performance is nearly independent of the choice of 
algorithm, since it is always within specification for an unfailed system. 
The relative merits of the different algorithms are then revealed by the 
system reliability. 

Still another definition might be: 

"Reliability is the probability that no missed 
alarms and less than three failures and false 
alarms combined (for the gyros, and the same 
for the accelerometers) occur during the 
mission. " 

With this definition, both reliability and performance depend upon the 
choice of algorithms. The logic of this choice is that false alarms, unless 
they lead to a missed alarm, are not likely to cause the large errors that 
missed alarms do. 

For the same algorithms, the first definition leads to higher relia- 
bility and worse performance than the third. Both the reliability and per- 
formance associated with the second definition depend upon the specifica- 
tion. The tighter the specification is, the higher the performance and 
the lower the reliability will be. 

The estimation of reliability requires the choice of one of these 
definitions that suits the mission and specification of the failure rate 
and failure magnitude probability density for each type of failure. The 
system error probability density for each failure magnitude must then 
be determined by simulation for each algorithm for the zero, one, and 
two failure cases, and all of this information must be combined in the 
appropriate manner to obtain reliability and performance estimates 
according to the chosen reliability definition. 

It was not possible to carry out such an elaborate program in the 
course of this study. Instead, the competing algorithms are all sub- 
jected to the same series of failures. The resulting system errors are 
then used as relative figures of merit to rank the algorithms. This 
technique avoids the detailed characterization of the failure modes 
peculiar to a particular system that would otherwise be required. 

For a manned mission such as Space Shuttle, performance is very 
important. That is, if two algorithms yield the same reliability, the one 
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that performs FDIC quickly before large velocity or orientation errors 
are accumulated is certainly to be preferred over a slower FDIC that does 
allow such errors to build up, in order to minimize the risk to the men 
on board. 

On a long -duration unmanned astronomical satellite, however, the 
errors building up during FDIC action are of relatively little importance, 
merely causing a temporary loss of observation time. Algorithms designed 
for such a mission may score somewhat poorly in performance compared 
to algorithms designed for a shuttle mission. However, it is felt that 
even for such algorithms it will be meaningful, or at least interesting, 
to look at their performance relative to the other algorithms. 

1. 3. 1 Reliability With Perfect FDIC 

If we assume that the FDIC algorithms are perfect, then the relia- 
bility is simply the probability that less than three gyro failures and less 
than three accelerometer failures occur during the mission. If p is the 
probability of a given instrument failing during a specified time interval, 
then the probability of m out of n instruments failing in the specified time 
interval is given by P(m, n) 


P(m, n) = 


n! 

m! (n-m) ! 


m, . ,n-m 

P (1-p) 


(1-D 


When a constant failure rate, X, and a time interval, At, are specified, 
the probability of a given instrument and its associated electronics fail- 
ing during the time interval is 


P = 1 


-\At 

- e 


( 1 - 2 ) 


The gyro or accelerometer reliability is given by 


2 

R = X p (™> 6 ) 

m = 0 


(1-3) 


for the dodecahedron configuration, since three failures cannot be isolated 
without external information. If R^ is the gyro reliability and R & is the 
accelerometer reliability, calculated according to Eqs. 1-1 through 1-3, 
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then the system reliability is 


(1-4) 


R = R R 
g a 

If there are any series elements (common power supplies, clocks, heaters, 
computer, etc. ) with a significant failure rate, the actual reliability will 
be much less than Eq. 1-4. 

1.3.2 Performance With Perfect FDIC 

The performance of instruments in the dodecahedron configuration 
can be calculated in a simple manner if it is assumed that a perfect FDIC 
combines the outputs of the unfailed instruments by least- squares, and 
if it is assumed that the errors of the unfailed instruments are statisti- 
cally independent with zero mean and unit variance. 

If the FDIC output errors are expressed in terms of the strapdown 
package coordinate system, they are correlated, and they differ for each 
different combination of failed instruments. Furthermore, generally the 
axis of greatest error will not lie along one of the coordinate axes, so 
that the greatest error will not be immediately apparent. Therefore, the 
covariance matrix for any given combination of failures is diagonalized 
by means of an orthogonal transformation to a new error coordinate 
system. The error coordinate system, in which the FDIC output errors 
are uncorrelated, can be related to the input axes of the failed and unfailed 
instruments so as to provide a physical picture of its orientation. The 
principal axes of the error ellipsoid lie along the axes of the error coor- 
dinate system, which are different for each different combination of 
failed instruments. 

For zero instruments failed, the error coordinate system axes may 
be taken in any direction (retaining orthogonality, of course). The error 
ellipsoid is a sphere. 

For one instrument failed, axis A is along the input axis of the 
failed instrument. Axes B and C may have arbitrary orientation about 
A. The error ellipsoid is a prolate spheroid with major axis along A. 

For two instruments failed, their input axes form two acute angles 
and two obtuse angles. Axis A bisects the acute angles. Axis B bisects 
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the obtuse angles. Axis C is perpendicular to both input axes. The error 
ellipsoid has its major axis along A and its minor axis along C. 

For three instruments failed, there are two equally probable possi- 
bilities. Either the three faces of the dodecahedron normal to the input 
axes of the unfailed instruments meet at a vertex of the dodecahedron 
(case A) or the three faces of the dodecahedron normal to the input axes 
of the failed instruments meet at a vertex of the dodecahedron (case B). 
Axis A passes through the center of the dodecahedron and the vertex. 

Axes B and C may have arbitrary orientation about A. For case A, the 
error ellipsoid is an oblate spheroid with minor axis along A. For case 
B, the error ellipsoid is a prolate spheroid with major axis along A. 

Average errors are computed by adding together all of the covari- 
ance matrices (without diagonalization) for a given number of failures 
and dividing by the number of matrices. The result turns out to be equal 
to the identity matrix multiplied by a scalar, the square root of which is 
called the average error. Thus, when averaged over all possible com- 
binations, the error ellipsoid for a given number of failures is a sphere. 
Table 1-1 presents the standard deviations of the errors in the error 
coordinate system and the average errors. 


Table 1-1. Errors With Perfect FDIC 


Number of 
Instruments 
F ailed 

Number of 
Combinations 

Error Coordinate 
System Axes 

Average 

Error 

A 

B 

C 

0 

1 

0. 707 

0. 707 

0. 707 

0. 707 

1 

6 

1. 000 

0. 707 

0. 707 

0. 816 

2 

15 

1. 345 

0. 831 

0. 707 

1. 000 

3 A 

10 

0. 727 

1. 345 

1. 345 

1. 176 

3 B 

10 

3. 078 

0. 831 

0. 831 

1. 902 

3 A and B 

20 




1. 581 
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1.4 THE ALGORITHMS 


Eight different algorithms have been selected for competitive 
evaluation. The algorithms rejected are all earlier or alternative efforts 


of authors 

whose algorithms have been 

selected. In roughly chronological 

order the selected algorithms are 
1. Adaptive 66 

(Ephgrave) 3, 7 

2. 

Fifteen Threshold 

(Evans and Wilcox) 8, 9 ’ 19 

3. 

Squared Error 

(Gilmore, McKern, and 

4. 

Bayesian Decision 

Oehrle) 11, 12,13,14 
(Gully) 15, 16 

5. 

Theory 

Maximum Likelihood 

(Wilcox) 17, 18 

6. 

Minimax 

(Potter and 

7. 

Adaptive 72 

Deckert) 19, 20, 21 * 22 ’ 23 > 24 ’ 25 > 26 
(Chien) 27,28 

8. 

Sequential 

29 

(Eckelkamp and Schiesser) 


The algorithms are described briefly here and in more detail in Para- 
graphs 3. 2. 1 through 3. 2. 8. However, the reader must go to the refer- 
ences for the theoretical backgrounds and derivations of the algorithms 
which, in many cases, are quite lengthy. 

3 7 

The Adaptive 66 Algorithm by Ephgrave ' uses a weighted least- 
squares estimator to perform failure detection, isolation, and correction. 
The estimator estimates the instrument package input vector (accelera- 
tion or angular velocity). The residuals are used to estimate the variances 
of the instrument errors. The variances are used to obtain the weights 
for the estimator in an iterative fashion. 

17 18 

The Maximum Likelihood Algorithm by Wilcox ’ is very simi- 
lar but uses six additional test signals whose expected values are the 
instrument variances to start the iteration. 

The Bayesian Decision Theory Algorithm by Gully * performs 
failure correction by means of a least-squares estimator using only the 
instruments classified as unfailed. The residuals are tested against 
thresholds to perform failure detection and isolation. 
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The Squared Error Algorithm by Gilmore, McKern, and 
Oehrle^’ ^ performs failure correction by means of a least- 

squares estimator like that of the Bayesian Decision Theory Algorithm. 

The sum of the squares of the residuals (total squared error) is tested 
against a threshold to perform failure detection. The squares of the 
residuals are tested against the total squared error to perform failure 
isolation. 

The Fifteen Threshold Algorithm by Evans and Wilcox^’ ^ ^ per- 
forms failure correction by means of a least-squares estimator like those 
of the preceding two algorithms. Fifteen test signals, each involving 
the outputs of a different set of four instruments, are compared to thresh- 
old to perform failure detection and isolation. 

The Minimax Algorithm by Potter and Deckert^’ ^ 

offers two different methods of failure correction. One is a least-squares 
estimator like those of the preceding three algorithms. The other is a 
Bounding Sphere Algorithm (described below). Failure detection and iso- 
lation are performed with the same fifteen test signals as the preceding 
algorithm but different logic. 

29 

The Sequential Algorithm by Eckelkamp and Schiesser uses a 
Kalman-Bucy filter to perform failure correction. The filter has nine 
states: three for the package input vector and six for the instrument 
errors. The measurements are the six instrument outputs. The K-B 
residuals and the instrument error states are compared against thresh- 
olds to perform failure detection and isolation. Failed instruments are 
eliminated from the measurement vector. 

27 28 

The Adaptive 72 Algorithm by Chien ' performs failure correc- 
tion by means of a least-squares estimator like those of the Bayesian 
Decision Theory, Squared Error, and Fifteen Threshold Algorithms. 
However, it also performs "identification" to decide whether a failure 
has a bias, a ramp, or a "variance" type of waveform. It then performs 
"recompensation" to estimate the value of a bias or ramp failure and 
"recertification" to decide if the recompensation effort was successful 
enough to allow the instrument to be reinstated using the new compensation. 
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It is reinstated simply by reverting to the least- squares solution that 
includes it. Failure detection and isolation use six test signals that 
are a subset of the fifteen test signals mentioned above for the first fail- 
ure and five test signals that are also a subset of the fifteen test signals 
for the second failure. The selection of five signals depends upon which 
instrument failed first so that all fifteen signals are used at one time 
or another. 

1.5 THE SIMULATION 

Figure 1-1 shows the overall system block diagram of the FDIC 
simulation program, FAILSIM. It is used to compare the different 
algorithms. 

• The trajectory generator provides an ideal trajectory that 
includes rotational and translational vibratory motions to 
exercise the algorithms. 

« The instrument configuration calculation provides compo- 
nents of angular velocity and acceleration along the input, 
output, and spin axes of each gyro and along the input, pen- 
dulous, and output axes of each accelerometer. 

• The unfailed errors calculation provides all of the errors 
experienced by the gyros and accelerometers in unfailed 
operation, except those produced by the sampling and 
quantization processes. 

• The additive failures calculation provides the errors caused 
by those failures that add to the output produced by the 
instruments in unfailed operation. (An example would be 

a bias drift shift which would leave the response to inputs 
unchanged. ) 

• The integration, sampling, and quantization calculations 
provide these functions, converting angular velocity and 
acceleration into increments of angle and increments of 
velocity. 

• The substitutional failures calculation provides the instru- 
ment outputs caused by those failures that substitute an 
incorrect output for the correct one (an example would be 
zero output). 
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© The nine failure detection, isolation, and correction (FDIC) 
algorithms each combine the six gyro outputs into an incre- 
mental angle vector and the six accelerometer outputs into 
an incremental velocity vector. FDIC algorithm 0 is a 
nominal algorithm supplied with information from the fail- 
ure simulator so that it always uses all of the unfailed 
instruments and none of the failed instruments. FDIC 
algorithm 0 provides a standard against which the candi- 
dates can be compared. 

• The nine identical strapdown algorithms convert incremen- 
tal angles and incremental velocities into actual orientation 
and actual inertial velocity. 

• The orientation and velocity errors are calculated using 
the ideal orientation and ideal inertial velocity from the 
trajectory generator. 

• Missed and false-alarm statistics, orientation errors, 
velocity errors, and computer time are calculated to pro- 
vide a basis for comparing the different algorithms. 

A computer session wherein the program is loaded and executed is 
referred to as a run. Each run may contain one or more cases, each 
case having its own input data. Each case may contain one or more 
Monte Carlo trials, all trials having the same input data but having 
different pseudorandom number sequences. 
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2. CONCLUSIONS 


Seven of the eight algorithms were simulated (all except the Adaptive 
72 Algorithm). The detailed results are given in Section 8. The con- 
clusions are based upon the gyro results only, since the accelerometer 
results were incomplete due to limited resources. Four of the algorithms, 
the Fifteen Threshold, Squared Error, Bayesian Decision Theory, and 
Minimax gave considerably better all-around performance than the other 
three. Adaptive 66, Maximum Likelihood, and Sequential. The first four 
have in common a definite logical structure that is absent from the Adap- 
tive 66 and Maximum Likelihood Algorithms and is not as clear-cut in 
the Sequential Algorithm. Thus we conclude that a definite logical struc- 
ture is to be desired. The structure should make explicit what instru- 
ments are classified as failed or unfailed, what hypotheses are being 
tested (or what decisions are being made), and what test signals and 
thresholds are to be used in testing each hypothesis (or making each 
decision). 

Of the three poorer algorithms, it seems clear that the Adaptive 66 
is the best and Maximum Likelihood the worst. The problem with the 
Maximum Likelihood Algorithm seems to be that subtracting off the residu- 
als from the inconsistency states is an inadequate method of preventing 
interaction between subsequent failures (Paragraph 3. 2. 5). The Sequen- 
tial Algorithm is also subject to interaction to a lesser extent (Para- 
graph 3. 1. 2.4). We conclude that interaction should be completely elim- 
inated as it is in the four better algorithms. (It is essentially eliminated 
in the Adaptive 66 Algorithm. ) By elimination of interaction we mean that 
the output of an instrument which has been classified as failed should not 
be used in making a decision about the classification of another 
instrument. 

Of the four better algorithms, none demonstrates a clear-cut super- 
iority. In the double simultaneous failure runs they rank as follows: 

Minimax 

Fifteen Threshold 

Squared Error 

Bayesian Decision Theory. 
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However, the Fifteen Threshold Algorithm is significantly poorer than the 
others in response to first failures. Apparently this result occurs because, 
unlike the other three algorithms, its thresholds for first and second fail- 
ures are identical. However, false alarms are more likely with one 
instrument already failed, requiring a higher threshold level for second 
failures than for first failures. Thus we conclude that independently 
adjustable thresholds for first and second failures (Squared Error, Baye- 
sian Decision Theory) or a properly preset ratio between the thresholds 
for first and second failures (Minimax) should be provided. 

The first and second failure plots show no clear-cut advantages for 
any of the remaining three algorithms, except possibly that the Bayesian 
Decision Theory Algorithm is not quite as good as the other two. If we 
let the double failure results govern, then the overall standings would be: 

Minimax 

Squared Error 

Bayesian Decision Theory 

Fifteen Threshold 

Adaptive 66 

Sequential 

Maximum Likelihood. 

However, differences between the top three algorithms are marginal and 
the ordering should not be taken too seriously. For example, the Mini- 
max Algorithm has relatively high computer sizing and timing require- 
ments. Weighting these factors heavily could move it down the list. 

We seem justified in making at least a tentative conclusion that 
indirect test signals are better than direct test signals (Paragraphs 3. 1. 2. 1 
and 3. 1.2.3), especially for double failures. This conclusion tends to 
substantiate the claims of Reference 24. 
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3. CANDIDATE FDIC ALGORITHMS 


The analytical background of the different algorithms is explored in 
Subsection 3. 1, while Subsection 3. 2 specifies the algorithms in detail. 

3. 1 ANALYTICAL COMPARISON OF ALGORITHMS 

The different algorithms are compared analytically in this section. 

3. 1. 1 General Considerations 

At least three types of information may be available for use by an 
FDIC algorithm: 

1. The outputs of the inertial instruments 

2. The internal states of the inertial instruments (such as 
gyro wheel speed) 

3. The outputs of electromagnetic radiation sensors (such 
as star trackers, doppler radars, etc. ). 

All of the algorithms use the first type of information. Only the 
Fifteen Threshold Algorithm uses the second type of information, although 
the others are capable of modification with varying degrees of facility to 
use it. In order to avoid undue complexity and keep the comparison fair, 
this capability will be deleted from the Fifteen Threshold Algorithm. 

None of the algorithms use the third type of information, but the Sequential 
Algorithm is formulated so as to be able to accept such information readily. 

There are at least two different concepts of the nature of the problem 
being solved when one designs an algorithm for this type of system. They 
lead to different types of algorithms. The first concept is that an instru- 
ment is either failed or unfailed. The algorithm must determine which 
instruments are failed and which are not. The subset of unfailed instru- 
ments is used for navigation while the failed instruments are ignored. 

The Fifteen Threshold, Squared Error, Bayesian Decision Theory, 
Minimax, Adaptive 72, and Sequential Algorithms are of this type. 

The second concept is that the important state of an instrument is 
the error of its output. The algorithm must determine the variance of 
this error and weight the instrument outputs accordingly. The Adaptive 
and Maximum Likelihood Algorithms are of this type. 
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The first problem concept may be considered as a special case of 
the second, in which the instrument error variance estimates can take on 
only the values one and infinity. 

The unfailed instruments are subject to errors during normal opera- 
tion. These errors can cause the algorithms to make mistakes in isolating 
the failed instruments. Both false alarms (unfailed instruments classified 
as failed) and missed alarms (failed instruments classified as unfailed) 
can occur. For the first problem concept, false and missed alarms appear 
as the use of an incorrect subset of instruments. For the second problem 
concept, false and missed alarms appear as the use of incorrect instrument 
weights. 

In the deterministic case, when the unfailed instruments have zero 
errors, it is possible to detect one, two, or three failures and to isolate 
one or two failures. However, when the errors resulting from the failures 
of two instruments have certain ratios, it is not possible to distinguish 
between that event and the failure of a different pair of instruments. (See 
Paragraph 3. 1. 3. ) In the deterministic case, this situation has probability 
zero, but when the errors of the unfailed instruments are nonzero, the 
probability of ambiguity becomes greater than zero. For this reason, the 
isolation of double simultaneous failures is difficult. Most of the algorithms 
seem to have been constructed under the assumption that double simultane- 
ous failures are very unlikely. Only the Minimax Algorithm provides a 
choice of implementations that depends upon this point. For consistency, 
only nonsimultaneous failures will be assumed for the Minimax Algorithm. 
(This decision does not preclude determination of the effects of double 
simultaneous failures on the algorithms during evaluation. ) 

Another question is how the different algorithms respond to "glitches" 
(temporary malfunctions), that is, whether or not an instrument will be 
reinstated as unfailed if "healing" is observed. The Fifteen Threshold and 
Minimax Algorithms do not permit reinstatement. The Adaptive 66, 

Squared Error, and Maximum Likelihood Algorithms do permit reinstate- 
ment. It seems from the documentation that the Bayesian Decision Theory 
Algorithm is not intended to permit reinstatement. If reinstatement is not 
permitted, then glitches are treated as failures. On the other hand, if 
glitches are infrequent in a particular system, prohibiting reinstatement 
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will preclude the erroneous reinstatement of failed instruments. Still 

another point is that reinstatement may permit the eventual correction of 

false alarms. A particular kind of glitch is the shift of bias drift rate 

sometimes experienced by gyros. The ability to recalibrate the drift and 

reinstate the instrument might be desirable. Such a capability, however, 

would introduce the possibility of the occurrence of unnecessary, erroneous 
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recalibrations. Gully discusses drift calibration but does not present an 
FDIC algorithm incorporating this feature. The Adaptive 72 Algorithm has 
the capability of recompensating and reinstating instruments with ramp 
errors as well as bias shifts. 

In the Adaptive 66 and Maximum Likelihood Algorithms, failure 
detection and isolation are obscured by a variance estimation scheme. In 
the Fifteen Threshold, Bayesian Decision Theory, Adaptive 72, and Sequen- 
tial Algorithms, failure detection occurs only when isolation occurs. In the 
Squared Error Algorithm, detection and isolation are separate processes. 
No action is taken when detection without isolation occurs. However, if an 
instrument is isolated as failed and subsequently the isolation test fails, 
the instrument is not reinstated unless the detection test fails also. The 
Minimax/ Bounding Sphere Algorithm has separate failure correction 
strategies for isolated and for detected but unisolated failures. 

The Fifteen Threshold, Squared Error, Bayesian Decision Theory, 
and Adaptive 72 Algorithms use a least-squares failure correction tech- 
nique. The Adaptive 66 and Maximum Likelihood Algorithms use a 
weighted least- squares failure correction technique. The Minimax 
Algorithm has a choice of either a least- squares failure correction tech- 
nique or the Bounding Sphere Technique, which minimizes the maximum 
possible estimation error. The Sequential Algorithm uses a Kalman-Bucy 
filter to perform failure correction. 

It is possible to construct six test signals, each predominantly 
dependent upon the error of a single instrument, but also dependent to a 
lesser extent upon the errors of the other instruments. Such "direct" test 
signals may be linear or quadratic and are used by the Adaptive 66, 

Squared Error, Bayesian Decision Theory, and Maximum Likelihood 
Algorithms. The linear direct test signals are proportional to the resi- 
duals of the least squares solution. 
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It is possible to construct 15 test signals, each depending upon the 
errors of a different group of 4 instruments, and completely independent 
of the errors of the other two instruments. Such "indirect" or "parity" 
test signals are used by the Fifteen Threshold, Minimax, and Adaptive 72 
Algorithms. 

The Fifteen Threshold, Squared Error, and Maximum Likelihood 
Algorithms all make use of prefilters to reduce the effects of quantization 
errors. To keep the comparison fair, the Adaptive 66, Bayesian Decision 
Theory, and Minimax Algorithms will be equipped with prefilters too. The 
Adaptive 72 Algorithm has a filter in its detection system and thus requires 
no prefilter. The Sequential filter is a modified Kalman-Bucy filter and 
also needs no prefilter. 

Since there are six instruments of a given type, six failed/unfailed 
decisions must be made or six variances must be determined to select or 
weight the instrument outputs properly. The six instrument outputs can 
provide six equations. However, these six equations contain nine unknowns: 
the six instrument errors and the three components of the angular velocity 
or acceleration inputs to the instrument package about which little a priori 
information is available. It is this deficiency of three equations that causes 
the fundamental difficulty in the design of failure detection algorithms 
based solely upon the instrument outputs and, by precluding the possibility 
of perfect failure detection, makes the problem interesting. 

Additional a priori information can be obtained by assuming that the 
magnitudes of the errors of four of the six instruments are less than 
quantities determined by their unfailed performance, or that the variances 
of four instruments must be equal to their unfailed value, while the vari- 
ances of the other two instruments may be greater than or equal to this 
value. The inequalities thus obtained are, however, insufficient to replace 
the missing equations. 

It is possible to solve for three "inconsistency states" which are 
linear combinations of the instrument outputs and therefore of the instru- 
ment errors and which are independent of the inputs to the instrument 
package and of each other. The vector comprising these three states is 
zero when all of the instrument errors are zero, regardless of the motion 
of the instrument package. This vector contains all of the failure detection 
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information available to the types of algorithm considered here. The 
various test signals used by the different algorithms are all functions of 
the inconsistency state vector. The only exception is the Sequential 
Algorithm which attempts to use a priori information about the vehicle 
motion. It is shown in Paragraph 3. 1.2.4 that this a priori information 
has little effect. 

Some of the algorithm documentation discusses running the algorithm 
at a slower rate than the rate at which information from the instruments 
is incorporated into the strapdown algorithm. For the purposes of this 
study, such an arrangement is considered unacceptable, because it permits 
a hard failure to wipe out the orientation or inertial velocity states com- 
pletely. Therefore all algorithms will run at the same sampling period, 
which is the strapdown algorithm minor cycle time. Thus they may be 
compared on an equal footing. (Of course, applications exist where a 
wipeout is acceptable. In general, they will require less computer 
capacity. ) 

3. 1. 2 Derivation of Equations 

The instrument outputs are: 

y = Ax + e (3- 1 ) 


where 

y = 6 -vector of instrument angular velocity or 
acceleration outputs 

A = 6 by 3 matrix of instrument input axis direction 

cosines with respect to axes fixed in the strapdown 
package 

x = 3-vector of instrument package angular velocity or 
acceleration, and 

e = 6-vector of instrument errors. 
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For all of the algorithms, we take 


c s 0 

c -s 0 

0 c s 

0 c -s 

s 0 c 

- s 0 c 


where 


c 


5 + vT 
10 


= cos a 


(3-2) 


(3-3) 


s = = sin a (3-4) 

^ : 

where a is half the angle whose tangent is two. The rows of A are unit 
vectors perpendicular to the faces of a dodecahedron, as seen in 
Figure 3-1. The columns of A are orthogonal 6 vectors of length \f2. It is 
possible to choose three more 6 vectors of length VZ which are orthogonal 
to the columns of A and to each other. Let them comprise the rows of a 
matrix C. One of the infinite number of possibilities is 


By definition 



(3-5) 


CA = 0 


(3-6) 




It is interesting to note that the ratio s/c is equal to the golden mean. 
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Figure 3-1. Dodecahedron Instrument Configuration 
Let v be a 3-vector of inconsistency states, where 

v = Cy (3-7) 

From Eqs. 3-i, 3-6, and 3-7 

v = C c (3-8) 

Thus the inconsistency states are independent of the IMU input angular 
velocity or acceleration and depend only on the instrument errors. All of 
the test signals used in the various algorithms can be expressed as linear 
or quadratic combinations of the components of v. 

Most of the algorithms have a prefilter which reduces the effect of 
high frequency noise, such as quantization error, on performance. The 
prefilter is a first-order filter such as 


f(i +1 )=k A f( i ) + kB v < i+1 > 


(3-9) 
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where 


A = exp (-T/t £ ) 

(3-10) 

: B = 1 " k A 

(3-11) 


with T being the sampling period and t £ the filter time constant. 

When the failures have been detected and isolated, it is necessary 
to correct for them. In the Adaptive and Maximum Likelihood Algorithms, 
correction is accomplished by obtaining the weighted least-squares 
estimate of x. 

x = B y (3-12) 

B = (a T Q" 1 a) _1 A T Q' 1 (3-13) 

where Q is the covariance matrix of the instrument errors and is assumed 
to be diagonal. 

In the Fifteen Threshold, Total Squared Error, and Bayesian 
Decision Theory Algorithms, Q takes on only the values 0 or 1, Thus, 
Eq. 3-13 becomes 

B=(a T a)“ 1 A T (3-14) 

where the rows of A corresponding to zero are set equal to zero. 

3. 1. 2. 1 Direct Test Signals 

Several of the algorithms make use of the residuals of the least- 
squares solution, as "direct" test signals, using them as estimates of the 
instrument errors. They are given by 

e =y - y (3-15) 

where 

y = A x (3-16) 

so that, from Eqs. 3-12, 3-15, and 3-16, 

* = (I - AB) y (3-17) 
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For the case of no failures, Q is the identity matrix and 


1 T 
B =■* A 


(3-18) 


One can see by direct calculation that 


1 TIT 

i-±aa t =ic T C 


(3-19) 


so that the c . are linear combinations of the v^: 

^ 1 T „ 1 _ T 

€ = 2 " C Cy = 2 C v 


(3-20) 


Thus the six residuals contain no more information than the three incon- 
sistency states. Prefiltered values of c can be obtained from 


- _ 1 r T , 
' =jC f 


(3-21) 


where f is given in Eq. 3-9. This formulation requires the prefiltering 
of only three quantities, rather than six. 
til 

When the k instrument has failed, the residuals are given. by 

til 

Eqs. 3-14 and 3-17 with the k row of A set equal to zero to give A(k). 

€ = {i - A(k) [a ( k) T A(k)] A(k) T }y (3-22) 

One can see, by direct calculation, that Eq. 3-22 is equivalent to 


6 =1 C(k) T C y C(k) T v 


(3-23) 


where 


c ij ,k) = c ij-l 1 ? 1 c <j c ik / / c i k 


(3-24) 


with the exception that the residual of the failed instrument, is given as 


€ k =y k 


(3-25) 
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by Eq. 3-22 and as 


=0 (3-26) 

k 

by Eq. 3-23. The latter is better, since the residual of the failed instru- 
ment has to be ignored, and Eq. 3-26 causes it to be ignored automatically. 


For the case of zero failures, the residuals may be calculated from 
Eqs. 3-1, 3-5, 3-6, and 3-20 
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(3-27) 


where 

d = — (3-28) 

v5 

The Bayesian Decision Theory Algorithm uses the residuals as test 
signals. 

3. 1. 2. 2 Variance Estimators 

The residuals, e ^ , may be considered to be estimates of the 
instrument error values. Their squares may be considered as estimates 
of the variances of the instrument errors. The squared residuals are 
used in the Squared Error Algorithm. The squared residuals are, how- 
ever, not unbiased estimates of the instrument error variances. From 
Eqs. 3-8 and 3-20 


A 

e 


1 T 

= I C Ce 


(3-29) 
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Let 


T 

D = C C 


then 


A 





e 


j 


and 


e . 

1 


6 6 

41 1 

j=l k=l 


D, ,D,. £ 
ij ik j k 


The covariance matrix of e is assumed to be diagonal, 

<Tk> = q i 4 


(3-30) 


(3-31) 


(3-32) 


(3-33) 


If we take the expectation of Eq„ 3-32 and introduce Eq. 3-33, we obtain 



1 _ 

4 



D 2 . q. 
U J 


(3-34) 


From Eqs. 3-5, 3-30, and 3-34, we obtain 



(3-35) 


and see that the estimates are biased. Unbiased estimates may be obtained 
by subtracting a portion of the total squared error from each squared 
error (residual). Let q be the unbiased estimate. The total squared error 
is 


6 

TSE = ^ e 2 . (3-36) 

j=l J 
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We now multiply the TSE by b and subtract it from the squared error. The 
result is normalized by a. 



To make q unbiased, b = 1/10, and to get the scale factor correct, a = 5. 


Thus 

l — / ■ * o 

6 

r — ' 



/v - ~2 1 \ /s2 

(3-40) 


V 

tl 

43 

M- 

(3-41) 

It is of 

some interest to express Eq. 3-40 in terms of the inconsistency 

states. 

Vj. From Eq. 3-20, we have 
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C i = 2 ^ C ji V j 

(3-42) 

and 

3 3 

5 i = * 2 Z C ji C W v j v k 

(3-43) 
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so that 


3 3 6 3 3 

s-iyTccvv >f 1 1 c c - - 

-i ' 4 pi k~i j 1 m > k 8 pi pi Pi 'J i _ki ' j -k 


t'l A A \ 
v ^--x-x; 


It can be shown from Eq. 3-5 that 


o 

> CL. C. . =2 6.. 

& jl ki J k 


so that 


3 3 


S li ( 5 ° k ‘ ‘ Z) 


v i T k 


Now if Q is a diagonal matrix 


Q u = 5 ij q i 


then Eq. 3-46 can be shown to be the solution of 

~ T T 

CQ C = v v 


which defines the test signals of the Maximum Likelihood Algorithm. 

a>Z A ^2 

Now let us compare e ^ with q^. For example, the value of t ^ is 

given by 


2 _ 1 T 
L " 4 e 
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*1 =€ 
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-28. 

The value of q^ 

is 
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(3-50) 


It is interesting to note that, if e k is nonzero and the other are 
all zero, then because of the zeros on the diagonal of Eq. 3-50 and the 

analogous equations for the other q, , only q is nonzero. This is not true 

, - - 2 
for or Cj . 


If we assume gaussian distributions, we can find the mean square 
error in the estimate of the variance. Assume that or^ is the variance of 
the k th instrument and cr 2 is the variance of the other five instruments; 
then it can be shown that 



Thus, both types of quadratic test signals have the same mean square 
error when considered as estimators of the variance of the instrument 
errors. Similar calculations can be done for the case of one failure. 

3. 1. 2. 3 Indirect Test Signals 

We now consider the "indirect" or "parity" test signals which, 
instead of trying to put in evidence the errors of a particular instrument. 
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show the error of a group of instruments that excludes one or two particu- 
lar instruments. Test signals that exclude one instrument are perfectly 
possible, but are not used by any of the candidates and will be ignored here. 
The Fifteen Threshold, Minimax, and Adaptive 72 Algorithms use the 
15 test signals, each of which excludes two instruments. Requiring that 
the signals be independent of the input angular velocity or translational 
acceleration determines them to within a constant factor. If u is the 
15-vector of test signals, then 


where a typical D matrix is 
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(3-53) 


(3-54) 
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Each component of u is a linear combination of the components of 
the inconsistency state vector v. Thus the matrix D can be factored into 
two matrices 

D = F C (3-55) 


and 


u=FCy=Fv (3-56) 

Prefiltering may be performed on the three components of v rather than 
the 15 components of u. Thus, we replace Eq. 3-56 by 

u = F f (3-57) 

where f is given by Eq. 3-9 and where 


0 

2 

0 

-s/c 

-c/ s 

1 

s/c 

-c/ s 

1 

-s/c 

c / s 

1 

s/c 

c/s 

1 

2 

0 

0 

-1 

s/ c 

c/ s 

1 

- s/c 

c/ s 

c/ s 

-1 

s/ c 

c/ s 

-1 

-s/c 

-1 

-s/c 

c/ s 

1 

s/c 

c/ s 

c/ s 

1 

-s/c 

c/ s 

1 

s/ c 

0 

0 

2 


31 



The Adaptive 72 Algorithm uses u^, u^, u^, u^q, u ll' an ^ u 15 ^ or 

detection and isolation of the first failure. For detection and isolation of 

the second failure when the first failure is instrument j, it uses the five 

u. which have D. , zero, 
i ij 

There is a close relationship between the indirect test signals and 

the least-squares residuals. If instruments j and k are assumed to be 

failed and a least- squares solution is performed excluding them, then all 

of the four residuals can be shown to be proportional to u^, where i is the 

index of the row of the matrix D for which D, , and D„ are both zero. Two 

ij ik 

of the residuals are equal to ±c/2 times u^, and two are equal to ±s/2 
times u^. 

3. 1. 2. 4 Sequential Algorithm 

The Sequential Algorithm is essentially a Kalman-Bucy filter, and 
thus appears rather different from the other algorithms. This section will 
use a coordinate transformation in 6 -space to demonstrate that the actual 
operation of the Sequential Algorithm is similar to the operation of the 
other algorithms. The detailed description is presented in Paragraph 3.2. 8. 

The linear, stochastic, multistage process to be filtered is 


x i+1 = $ Xj + u. (3-59) 

where 

<x o >= 0 (3-60) 

<^u^)>= 0 (3-61) 

<x x M (3-62) 

\o o / o 

<0*4 a l T )= Q 6 ij (3 ' 63) 

<u. x o T V 0 (3-64) 

Measurements y^ are linearly related to the state by 

y^ = HjX^ + w i (3-65) 
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where 


(3-66) 

(3-67) 

(3-68) 

(3-69) 


The filter is given by the following. For resets: 

=x. + K t (y. - (3-70) 

K i = M i H i T ( H i M i H i T + R ) (3-71) 

P 4 = (I - K i H i jM i (3-72) 

For propagation between resets: 

x i+ i = $ Xj (3-73) 

M, , , = $ P.$ T + Q (3-74) 

i+ 1 l 


The state vector, x, is a 9-vector whose first three elements are the 
package input 3-vector of acceleration or angular velocity. The last six 
elements are the instrument error 6 -vector. The measurement vector, y, 
is a 6-vector whose elements are the instrument outputs. Thus the 
measurement matrix is, from Eq. 3-1 

H =[a;I 6 ] (3-75) 

where A is given by Eq. 3-2. 
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The state transition matrix is 


$ = 




' L D 


The covariance matrices are 


M = 
o 


m l X 3 


Q = 




m 2 r 6 


0 


R = 


' [^ 6 ] 


(3-76) 


(3-77) 


(3-78) 

(3-79) 


Because of the H matrix, all of the nine states are coupled together, 
making insight into the filter operation rather difficult. To remove this 
difficulty, we transform the six instrument error states and the six 
measurements (instrument outputs) into a new coordinate system by means 
of the matrix D. 



(3-80) 


(3-81) 


y 1 = D y 


(3-82) 


where G is given by Eq. 3-5. It can be shown that the inverse of D is 


D 


-1 



(3-83) 
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If we carry out the transformation, we find that the new filter is given by 
Eqs. 3-70 to 3-74 with unprimed quantities replaced by primed quantities 
and with 


$ ' = $ 

M» = M 
o o 

Q' =Q 
R 1 = R 

However, the new measurement matrix becomes 



(3-84) 

(3-85) 

(3-86) 

(3-87) 


(3-88) 


Now we see that in all matrices ($', M' o , Q', R', and H'), the first six 
states are decoupled from the last three. Furthermore, because of the 
identity submatrices, there are actually six uncoupled sets of states com- 
prising three identical two-state filters (states 1 and 4, 2 and 5, 3 and 6) 
and three identical one-state filters (states 7, 8, and 9). (Note that this 
conclusion will no longer hold after failure detection causes the H matrix 
to be modified by deletion of a row. ) 

For the two- state filters, the Kalman gains can easily be found to be 


m. 


K« = K' =K' 


11 22 33 m^ + m 2 + r 


m„ 


K' = K' = K' 


41 52 63 mj + m 2 + r 


while, for the one- state filters, the Kalman gains are 


m. 


K» = K' = K' 


74 85 96 m 2 + r 


(3-89) 

(3-90) 


(3-91) 
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Now, since represents the variance of the vehicle motions while 
represents the variance of the instrument errors and r represents the 
variance of the measurement errors (quantization), we have that 

> > m^ (3-92 ) 


m^ > > r 


(3-93) 


Therefore K'^, K',^, anc ^ ^*63 are near ty zero an< ^ states 4, 5, and 6 of 
our transformed filter remain essentially zero. States 7, 8, and 9 are 
simply the result of three first-order Kalman filters acting upon y^, y^, 
and y^. However, from Eqs. 3-80 and 3-82 we see that 

u) 

y 5/ = r Cy= T V . (3-94) 

y 6 ) 


so that the inputs to the three first-order filters are proportional to the 
familiar three uncertainty states. Thus x'^, x'g, and x'^ are similar to 
the filtered uncertainty states f, except for the time -varying nature of the 
Kalman filter. Furthermore, from Eqs. 3-81 and 3-83 we see that 



(3-95) 


Since x'^, x'^, an d x, £ are essentially zero, x^ - x^ are analogous to 
filtered versions of the direct test signals or residuals discussed previ- 
ously. The smallness of x'^, x'j., and x'g comes about because the a 
priori information available about the vehicle motions is insufficient to 
contribute significantly to the solution of the failure detection problem. 
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Of course, if other types of measurements were taken (i. e. , star tracker) 
the situation would be completely different. It is in the latter case that 
the Sequential Algorithm should prove most effective. Unfortunately, it 
is beyond the scope of this study. 

One aspect of the Sequential Algorithm in which it differs from most 
of the other algorithms should be mentioned. When a failure is detected, 
failure correction is accomplished by deleting a row from the H matrix 
and proceeding as before. However, no changes are made to the existing 
values of the instrument error states. Eq. 3-27 shows that, when one 
instrument exhibits an error, it not only contributes to its appropriate 
residual, but contributes an amount equal to \/\fS (44.7%) of the error to 
the other residuals. When a failure is detected, the other five instrument 
error states will therefore include such errors of up to 44. 7% of the 
threshold level, which will subsequently gradually decay with time. In 
most of the other algorithms, the FDIC reorganization after a failure is 
such that no part of the output of a failed instrument contributes to the test 
signals used to find the second failure. 

Since the Sequential Algorithm has so many adjustable parameters, 
a relationship has been derived between its Kalman-Bucy filter operating 
in steady state and the constant gain prefilters used by many of the other 
algorithms, in order to gain some insight into the meanings of the param- 
eters. This derivation is presented in Appendix B. 

3. 1. 3 Multiple Simultaneous Failures 

The possibility, mentioned at the beginning of this section, that two 
failures can be indistinguishable from two other failures in the determi- 
nistic case can be seen by considering the three sets of errors 

s T = (c 0 0 0 s 0) 

e T = (0 -c 0 0 0 s) (3-96) 

e T = (0 0 -c -s 0 0) 

For each of these different cases of two instruments failed and four 

instruments unfailed, multiplication of the error vector by the matrix C 

2 2 

of Eq. 3-5 gives v = (0, -c , s ). Since we have shown that all of the 
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test signals of aLl of the different algorithms are functions of v, the three 
cases are indistinguishable to all of the FDIC algorithms. Reference 24 
gives an interesting discussion of simultaneous failure detection and 
isolation. 

3. 1. 4 Linearity Versus Nonlinearity 

Occasionally it has been suggested that the instrument biases should 
be estimated by some kind of linear filter and then used to adjust the bias 
compensation of the instruments, thus avoiding the necessity of switching 
out the instruments or weighting their outputs. This procedure is linear. 

It is therefore doomed to failure. The least- squares solution already gives 
the optimal linear estimate of acceleration or angular velocity. For 
example, suppose we modify the instrument outputs by subtracting the 
estimates of the instrument errors (residuals) from them 

y 1 = y - € (3-97) 

where c is given by Eq. 3-15. Then an improved estimate might be 
obtained by modifying Eq. 3-12 to 

x' = B y' (3-98) 


From Eqs, 3-97, 3-98, and 3-15 


/s 

= B y 


From Eqs. 3-99 and 3-16 


B A x 


(3-99) 


(3-100) 


From Eq. 3-13 we see that B A is the identity matrix so that 

x« = x (3-101) 

and no improvement occurs. This behavior is a well-known feature of 
optimal linear solutions, intimately related to their optimality. The same 
situation arises even if a Kalman-Bucy filter is used to estimate the 
instrument errors, because the estimate of the package inputs is already 
the optimal estimate. 
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All of the FDIC algorithms studied here are nonlinear. The "linear 
parts" of all of them have been shown above to be essentially identical, 
all test signals depending upon the inconsistency vector v. It is the man- 
ner in which the necessary nonlinearity is introduced that distinguishes 
each algorithm from the others. 

From Eq. 3-19 we can see that the instrument error vector can be 
broken up into two parts by projection matrices 

IT IT 
e AA X € + j C ± C € 

These two vectors are orthogonal to each other in 6-space 

(i A A T .) T (|c T cA =i. T AA T C T C. 

Transposing Eq. 3-6 gives 

A T C T = 0 (3-104) 

So that Eq. 3-103 must equal zero. If we write Eq. 3-102 as 

€ = c + e (3-105) 

we see that e contributes to the estimation error 3-vector but not to the 
inconsistency 3-vector, and that e contributes to the inconsistency vector 
but not to the estimation error vector. « is accessible to the FDIC 
designer, but ? is not since the actual package motion is not known. The 
FDIC problem is to use knowledge of e to find e without knowledge of e. 
Fortunately, perfect knowledge of « is not required. 

3. 1. 5 V Space Logic 

The 3-vector v is defined in some 3- space which we may call v- space. 
This space is not the same as the ordinary 3-space in which the instrument 
input axis vectors and the package input vector x are defined. In the preced 
ing section, we see that v-space and ordinary 3- space may be considered 
as having orthogonal bases in a 6-space. 


(3-102) 


(3-103) 



It is interesting to think of the different FDIC algorithms as struc- 
tures in v-space. Since all of the test signals are functions of the v^, the 
comparison of a test signal to a threshold appears as the determination of 
whether the end of the v-vector lies on one side or the other of a surface 
in v-space. For example, total squared error equal to a threshold level 
generates a sphere centered on the origin of v-space, and a parity signal 
equal to a threshold level generates a plane. 

Thus the failure isolation algorithm synthesis problem can be looked 
upon as setting up the optimal boundary surfaces in v-space, associating 
decisions with the crossing of various boundaries, and changing the 
boundaries optimally after each decision. 

In an algorithm with a prefilter, the filtered value of v, f, would be 
used in place of v. 

Up to the present time, no one has performed FDIC synthesis directly 
in v-space. Perhaps the concept will never be of more than academic 
interest. 

3. 1. 6 Filtering 

The quantization error experienced with the typical strapdown 
package calls for some type of a filtering process to reduce its effect. 

For example, a 5-arc sec quantum size and a 10-msec sampling period lead 
to an instantaneous gyro drift rate with a standard deviation of 204 deg/hr. 

The Fifteen Threshold and Maximum Likelihood Algorithms make 
use of a simple first-order low-pass prefilter. First let us consider the 
effect of step failures. If we approximate the sampled-data filter by a 
continuous filter, its transfer function is simply 


G(s ) = 


1 

1 + T^S 


(3-106) 


Where is the filter time constant. The filter response to a step input 
of magnitude a is 


u f (t) = a 



(3-107) 
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and the time required to reach a threshold level, b, is 

t = t In (3-108) 

The error caused by a step (angular velocity or acceleration) of magnitude 
a is approximately 

e = at (3-109) 

(angle or velocity). 

Therefore 


t 


= a t. 


In 


a 

a-b 


If we normalize Eq. 3-110 on the threshold level b. 



(3-110) 


(3-111) 


Figure 3-2 shows this relationship. The larger the failure is, the less 
error it contributes to the system before it is switched out. If the failure 
is less than the threshold (a/b < 1), then the error increases indefinitely 
until the end of the mission. 


Now let us consider the effect of quantization errors. The transfer 
function of the sampled-data filter (of which Eq. 3-106 is an approximation) 
is 



(3-112) 


The angular velocity output of the gyro (or the acceleration output of 
the accelerometer) is given by continuous integration, sampling, quantiza- 
tion, and numerical differentiation as shown in Figure 3-3. 
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Figure 3-2. Low-Pass Prefilter Error 
Versus Failure Magnitude 
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Figure 3-3. Gyro or Accelerometer Model 


Now let us approximate the quantization error, e*, by a random 
sequence such that 



(3-113) 

(3-114) 


where the standard deviation is given by 


2 

<r 



(3-115) 
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where Q is the quantum size. (Eq. 3-115 assumes that e^ is uniformly- 
distributed over ±Q/2. ) From Eq. 3-112 and Figure 3-3, we see that the 
filter output caused by quantization is 



By long division, Eq. 3-116 becomes 



(3-118) 


From Eqs. 3-114, 3-115, and 3-118, 
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For T f — 0 


Q 

W f sfl T 


the unfiltered, value; while for ® 




(3-121) 


(3-122) 


the error becomes independent of the sampling period. From Eq. 3-120 
we may write 




(3-123) 


This function is plotted in Figure 3-4. 



NORMALIZED DRIFT RATE 


Figure 3-4. Normalized Drift Rate Versus Normalized Time Constant 
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The Squared Error Algorithm makes use of a different type of 
prefilter. The instrument outputs are accumulated. Every 60 sec the 
accumulators are purged of all outputs older than 120 sec. If we approxi- 
mate the sampled-data accumulator by a continuous one, the prefilter looks 

like Figure 3-5, where T is the purge period of 60 sec. (The gain of 

P 

1/3 T normalizes the prefilter so that its steady- state peak gain is unity, 

P 

to correspond with Eq. 3-107. This gain factor does not actually appear 
in the Squared Error prefilter and is used here only to facilitate compari- 
son with the low-pass prefilter described previously. ) 

The response of this prefilter depends upon the phasing of the failure 
with respect to the purges. The response to a step failure of magnitude a 
is shown in Figure 3-6 for four different phases. By averaging over all 
values of relative phase, we can find the mean time to exceed a threshold 
level, b. 

0 < a/b < 1 

1 < a/b < 3/2 (3-124) 

3/2 < a/b < ® 

0 < a/b < 1 

1 < a/b < 3/2 (3-125) 

3/2 < a/b < co 


t = 


3— + 2 
a 


IN 


3 — T 
a p 


Following Eq. 3-109, and normalizing. 


3bTp 


[!!--!!] 


This function is plotted in Figure 3-7. 



Figure 3-5. Squared Error Algorithm Prefilter 


45 








Figure 3-6. Prefilter 
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Figure 3-7. Squared Error Algorithm Prefilter 
Error Versus Failure Magnitude 

A comparison of Figure 3-7 with Figure 3-2 shows a more desirable 
characteristic for the Squared Error Algorithm prefilter. The error is 
always less than that of the simple first-order prefilter. 

Now let us consider the effect of quantization errors. Referring to 
Figure 3-5, the analog integrator is actually the pulse accumulator with 
Z transform Tz/(z-l). This transform cancels the one: in the right-hand 
box of Figure 3-3 so that 


' TT [ e(t) - e <‘- 2T P >] 


Using the same quantization error model as before, we see that 

<r ri- (3-12-; 

“f Jl 3T p 

Since we see from the normalization procedure that and 3T^ are 
approximately equivalent, Eq. 3-127 shows an error due to quantization 
about y/z times larger than Eq. 3-122. 
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These results do not provide grounds for a clear choice between the 
two prefilters, since the low-pass filter has the better response to quanti- 
zation error, and the Squared Error Algorithm prefilter has the better 
response to step failure. Any of these algorithms could use either 
prefilter. To be able to determine the relative merits of the algorithms 
themselves, the same prefilter will be used for all three. In addition, 
the Adaptive 66, Bayesian Decision Theory, and Minimax Algorithms, 
which contain no filtering at all, will be provided -with the same prefilter. 
The first-order, low-pass type of prefilter is selected for this purpose 
because it is simpler. After a single algorithm is selected, a prefilter 
tradeoff study can be made in the future. 

The Sequential Algorithm is essentially a modified Kalman-Bucy 
filter and therefore requires no prefilter. It is discussed in more detail 
in Paragraph 3. 1. 2. 4 and Appendix B. 

The Adaptive 72 Algorithm is considerably more elegant than any of 
the above in its filtering techniques. Originally it used a Kalman-Bucy 
filter to whiten the unfailed instrument errors prior to the detection system. 
Later, when the algorithm was implemented on SIRU, it was felt that the 
errors were close enough to white so as to make this filter unnecessary. 

The detection system uses a filter which can be represented as 
shown in Figure 3-8 if the digital accumulator is represented by an analog 
integrator. This detection system is for positive polarity failures; a 
similar system is used for negative polarity failures. 



Figure 3-8. Detection System for Positive Polarity 
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The loop about the integrator resets it to zero at every sampling 
instant if the integrator output is less than zero. The constant b sets a 
level of angular velocity error or acceleration error below which no fail- 
ures are detected. The constant c determines how quickly a failed instru- 
ment is switched out. The output f is 0 for no failure and 1 for a failure. 
For a step failure of magnitude a, the time required to reach the threshold 
level c is 


Following Eq. 


! co o 


0 < a/b < 1 
< a/b < as 


3-109, and normalizing. 


b(c/b) 


CO 

a/b 

a/b - 1 


0 < a/b £ 1 

1 < a/b < ® 


(3-128) 


(3-129) 


This function is plotted in Figure 3-9. The error lies above that shown 
in Figures 3-2 and 3-7. 


b(c/b) 



Figure 3-9. Adaptive 72 Algorithm Filter Error 
Versus Failure Magnitude 
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By analogy with the Squared Error Algorithm prefilter, we see that 
the equivalent quantization error response is 


Q 1 

T “f "7b =7b 


(3-130) 


From the normalization process c/b is approximately equivalent to 3T^, 
so that the quantization error response of the two systems is comparable. 
It should be kept in mind that the above comparison does not take into 
account the nonlinear characteristic of the Adaptive 72 filter, which 
should have an important effect on the detection of signals in the presence 
of noise. 

3.2 DETAILED ALGORITHM DEFINITIONS 


The algorithms, as they appear in FAILSIM, are specified in detail 
in this section. However, the reader is referred to the original docu- 
mentation of the algorithms for their explanations, derivations, or ration- 
ales, some of which are quite lengthy. Each algorithm seems to have 
been developed with a somewhat different set of assumptions about the 
nature of the problem being solved. 

3. 2. 1 Adaptive 66 Algorithm 

The Adaptive 66 Algorithm is described in References 3 and 7. It 
uses a weighted least-squares estimator to perform failure detection, 
isolation, and correction. The estimator estimates the instrument pack- 
age input vector (acceleration or angular velocity). The residuals are 
used to estimate the variances of the instrument errors. The variances 
are used to obtain the weights of the estimator. New estimates are then 
obtained in an iterative fashion until the variances cease changing. 

In the original algorithm, the equation for the variance in terms of 
the squared residuals is given by 

Q ii = *i + *o ' 1 = 6 (3-131) 

2 

where is the nominal or minimum value. The computation time can be 
reduced by several orders of magnitude if Eq. 3-131 is replaced by 
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since the Q_ in Eq. 3-132 will not change for small perturbations of the 
e. caused by unfailed instrument errors. If the do not change, it is 
unnecessary to recompute the correction matrix B. The change from 
Eq. 3-131 to Eq. 3-132 is necessary to prevent this algorithm from con- 
suming an inordinate amount of simulation computer time. In addition to 
this change, a prefilter has been added which filters the instrument out- 
puts prior to their use in the algorithm. Figure 3-10 shows the algorithm 

flow diagram. The Adaptive 66 Algorithm has two adjustable parameters, 

2 

the prefilter time constant t , and the a priori error variance e . 


Initialization 



Figure 3-10. Adaptive 66 Algorithm Flow Diagram 
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The flow diagram description follows: 

(1) The filter 6 vector_f, the previous value of the instrument 
covariance matrix Q, and the correction matrix B are 
initialized before the trial. 

(2) Each minor cycle the filter outputs are calculated. The 
iteration counter is zeroed. 

(3) The iteration count is increased by one. 

(4) If 100 iterations have occurred, the iteration loop is 
te rminated. 

(5) The residuals are calculated from the filter outputs, and 
a new covariance matrix is obtained. 

(6) If the changes in all of the diagonal elements of the 
covariance matrix from their previous values are less 
than the factor CHIT times their previous values, the 
iteration loop is terminated. 

(7) The previous value of the covariance matrix is set equal 
to the new value. The correction matrix is updated. A 
new iteration is performed. 

(8) The package input is estimated from the unfiltered 
instrument outputs. 

3. 2. 2 Fifteen Threshold Algorithm 

The Fifteen Threshold Algorithm is described in References 8, 9, 
and 10. It makes use of the fifteen indirect test signals described in 
Paragraph 3. 1. 2. 1. Each element of the test signal vector u is compared 
to a threshold. A vector w is used to store the results of this compari- 
son. If the threshold is exceeded by |uj, then w^ is set equal to 1. If 
|uj later becomes less than the threshold, then w^ is not reset to 0. Thus 
there is no provision for reinstatement of "healed" instruments. 

An instrument is considered failed only if all of the w^ to which it 
contributes a signal are equal to 1. Eqs. 3-12 and 3-14 are used for 
failure correction to save time in the simulation, although the original 
algorithm used Gauss- Jordan reduction so that, in a real-time system, 
a change to the failed /unfailed status could be implemented within a sin- 
gle minor cycle. Figure 3-11 shows the algorithm flow diagram. 
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Figure 3-11. Fifteen Threshold 
Algorithm Flow 
Diagram 
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The variables S and S' are 6-bit binary numbers whose bits 
correspond to the unfailed (0) or failed (1) states of the instruments. S 
is the instrument state. S' is set to all ones at the beginning of each 
cycle. Whenever a test signal is found to be less than the threshold, the 
1 bits corresponding to the instruments whose outputs contribute to the 
signal are masked out by the appropriate mask: 


1 

1 

2 

3 

4 

5 

6 

7 

8 
9 

10 

11 

12 

13 

14 

15 


AMASK. 

i 

0 0 0 0 1 1 

0 0 0 1 0 1 

0 0 0 1 1 0 

0 0 1 0 0 1 

0 0 10 10 
0 0 110 0 
0 1 0 0 0 1 

0 10 0 10 
0 10 10 0 
0 110 0 0 
1 0 0 0 0 1 

1 0 0 0 1 0 

1 0 0 1 0 0 

10 10 0 0 

1 1 0 0 0 0 


Those instruments whose bits are not masked out are failed. The state 
of a particular instrument is determined by masking out the states of the 
other instruments with the appropriate mask: 


OMASK. 

x 


1 

2 

3 

4 

5 

6 


1 0 0 0 0 0 
0 1 0 0 0 0 
0 0 1 0 0 0 
0 0 0 1 0 0 
0 0 0 0 1 0 
0 0 0 0 0 1 
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The test signals are prefiltered before use in the algorithm. The Fifteen 
Threshold Algorithm has two adjustable parameters, the prefilter time 
constant T r and the threshold level TH. 

x 

The flow diagram description follows; 

(1) The filter 3 vector f, the instrument state S, the threshold 
state 15 vector w, and the correction matrix B are initi- 
alized before the trial. 

(2) Each minor cycle the intermediate instrument states S' 
are set to all ones, and the inconsistency state vector v 
is calculated and filtered. 

th 

(3) The 15 threshold states are examined. If the i state is 
equal to 1, it is ignored. 

th 

(4) If it is equal to 0, the i parity signal is calculated. 

th 

(5) The i parity signal is compared to the threshold. 

(6) If the threshold is not exceeded, all of the bits correspond- 
ing to instruments included in the test signal are masked 
out to 0. 

(7) If the threshold is exceeded, the threshold flag is set 
equal to 1. 

(8) If S = S', no change in failure status has occurred, and 
the correction matrix calculation is bypassed. 

(9) If S' is all ones, more than two failures have been detected 
and isolation is not possible. S is not changed, and the 
correction matrix calculation is bypassed. Since it cannot 
resolve the situation, the algorithm continues to use the 
old correction matrix. 

(10) S is set equal to S', and the new correction matrix B is 
calculated. The new B matrix is calculated from Eq. 3-14 
with the rows of A corresponding to the failed instruments 
taken as zero. 

T 

(11) The matrix A A is zeroed. 

(12) The k instrument is tested for failure. 

tli T 

(13) If the k c instrument is qnfailed, its contribution to A A 

is calculated and added in. 

(14) (A^ A) ^ is calculated. 
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fh 

(15) The j instrument is tested for failure. 

(16) If the j^ 1 instrument is failed, the j^ 1 row of B is set 
equal to zero. 

x.-J_ 

(17) If not, the j row of B is calculated. 

(18) The package input is estimated from the unfiltered 
instrument outputs. 

3. 2. 3 Squared Error Algorithm 

The Squared Error Algorithm is described in References 11, 12, 13, 
and 14. The algorithm actually used in FAILSIM differs in several ways 
from these descriptions. The differences and the reasons for them will 
now be discussed. 

The prefilter used with the Squared Error Algorithm in SIRU con- 
sists of an accumulator into which the instrument outputs are summed. 
Every 60 sec, the accumulator is reset so that it contains only the last 
120 sec worth of accumulated output. Thus the data accumulation period 
fluctuates from 120 to 180 sec. The prefilter used in FAILSIM is a sim- 
ple first-order, low-pass filter, of the kind used for most of the other 
algorithms. The reason for this change is that the prefilter is not an 
integral part of the different algorithms. If the SIRU prefilter should be 
superior to the low-pass prefilter, all of the algorithms can make use of 
it. The purpose of FAILSIM is to compare the basic algorithms with each 
other. This can best be done if the prefilters are the same. After a 
single algorithm is selected, a prefilter tradeoff study can be made in 
the future. 

The Squared Error Algorithm automatically raises the failure 
detection thresholds for the gyros when the angular velocity becomes 
large, in order to reduce the probability of false alarms caused by instru- 
ment dynamic errors, scale factor errors, and misalignments. In SIRU, 
the signal used for this purpose is the square of the sum of the absolute 
values of the three components of angular velocity. This signal is not 
a scalar; that is, its value depends upon the choice of coordinate system. 
Since the choice of coordinate system is arbitrary, the signal is charac- 
terized by an undesirable quality of arbitrariness. It is more usual to 
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use invariant or covariant expressions, i. e. , scalars, vectors, and 
tensors. Therefore the squared magnitude of the angular velocity vector 
will be used for this purpose in FAILSIM, since it is a scalar. 

If a spike of angular velocity occurs and causes a spike of gyro 
errors, the threshold would be raised as described above during the spike 
and then lowered. The filtered inconsistency states, however, would 
exhibit a decaying exponential at the filter time constant, which could 
cause the now lowered threshold to be exceeded. (It would seem that the 
SIRU prefilter would also experience a form of this problem.) Therefore, 
instead of raising the squared error threshold, the FAILSIM algorithm 
lowers the gain at the input to the prefilter. Since the gain acts upon the 
unsquared signals, it is taken as the square root of the reciprocal of the 
gain applied to the squared error threshold in SIRU. 

The logic in the SIRU algorithm provides that, when one failure has 
been detected, the algorithm will alternate between testing for a second 
failure and testing for healing of the first failure, so that each of these 
two tests occur only every other minor cycle. This procedure is required 
because of the limited time available in the particular computer used for 
SIRU. However, since all of the other algorithms are allowed all the time 
they need, it seems more equitable to modify the logic so that both tests 
take place each minor cycle, and this action is taken* 

The SIRU algorithm provides checks for overflow of the error quan- 
tities. Since overflow is not a problem for the CDC 6500, these checks 
have been eliminated. 

The SIRU algorithm detects the third failure. However, it does 
nothing with the information other than output it. Therefore, this capa- 
bility is deleted to conserve computing time. 

The squared-error calculations have been modified to make the 
coding simpler. However, the results are numerically identical to those 
of the original equations. When zero instruments are failed, the algorithm 
makes use of the residuals as estimates of the instrument error. If pre- 
filtering is performed, we have from Eq. 3-21 

€=|c T f (3-133) 
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(3-134) 


The total squared error for zero failures is defined as 

TSE = € T € 

From Eqs. 3-5, 3-133, and 3-134 we have 

TSE q =|f T f (3-135) 

th. 

When the k instrument has failed, the residuals are calculated as in 
Eq. 3-23 

€ (k) = j C(k) T f (3-136) 

The total squared error for the remaining instruments is 

TSE k = € "(k) T e(k) (3-137) 

Note that because of Eq. 3-26, the residual of the failed instrument is 
automatically dropped from TSE^. 

The total squared error is compared to a threshold. If the thresh- 
old is exceeded, then the ratios of the individual squared errors (squared 
residuals) to the total squared error are tested against another threshold. 

If that threshold is exceeded by one of the ratios, the corresponding instru- 
ment is considered failed. The same procedure is used for the second 
failure. A flag LHEAL controls the healing process. 

LHEAL = 0 healing is not allowed 
LHEAL = 1 healing is allowed 

Figure 3-12 shows the algorithm flow diagram. If the flag LHEAL 
is zero, instruments are not reinstated even though healing may be 
observed. The Squared Error Algorithm has six adjustable parameters, 
the prefilter time constant t^, the threshold levels K^, K^, K^, and K^, 
and the gyro dynamics compensation gain 
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The flow diagram description follows: 

(1) The filter 3 vector f, the instrument state S, the failed 
status indicators I and J, and the correction matrix B 
are initialized before the trial. 

(2) Each minor cycle the filter gain is calculated as a func- 
tion of the previous angular velocity estimate (gyros only). 

(3) The inconsistency state vector v is calculated and filtered, 
the total squared error is found, and the present values 

of S and I are stored. 

(4) The total squared error is tested against threshold Kl. 

(5) If the threshold is exceeded, the residuals are calculated. 

(6) The ratios of the squared residuals to the total squared 
error are tested against threshold K2. If the threshold is 
exceeded, the loop terminates and control passes to Step 7. 
If no ratio exceeds the threshold, control passes to Step 16. 

(7) K is the number of the failed instrument just detected. I 
is the number of the "first" failed instrument, and J is 
the number of the "second" failed instrument. If I already 
equals K, no further action is required, and control passes 
to Step 10. 

(8) K is tested against J. If it equals J, then the values of I 
and J are exchanged in Step 11. 

(9) If not, I is tested against 0. If I equals 0, no failure was 
previously detected, and I is set equal to K in Step 13. 

(10) If not, J is tested against 0. If J equals 0, only one 
failure was previously detected. 

(11) Its instrument number is stored in J, and I is set equal 
to K. 

(12) If J is not equal to zero, then there were two previous 
failures, neither one of instrument K. If healing is not 
possible, no further action is taken. 

(13) If healing is possible, I is set equal to K and J is zeroed. 

(14) If healing is not possible, total squared error less than 
the threshold has no effect. 

(15) If healing is possible, any previous failed instruments are 
reinstated when the total squared error is less than the 
threshold. 
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(16) If I equals zero, no instruments are failed. 

(17) The instrument states are zeroed. 

(18) If I is not zero, the 1^ instrument state is set equal to 1. 

(19) If I is unchanged Step 20 is bypassed. 

(20) The matrix for calculating the residuals from the filter 
vector when instrument I is failed is calculated. 

(21) The residuals and total squared error (excluding instru- 
ment I) are calculated. 

(22) The total squared error is tested against threshold K3. 

(23) The ratios of the squared residuals to the total squared 
error are tested against threshold K4. If the threshold 
is exceeded, the loop terminates, and control passes to 
Step 25. If no ratio exceeds the threshold, control passes 
to Step 28. 

(24) If the total squared error is less than the threshold, K is 
set equal to 0. 

(25) If J equals 0, there has been no "second failure. " 

(26) If J is not equal to zero, then the healing flag is tested. 

(27) If J equals 0 or healing is permitted, J is set equal to K. 

(28) If J equals 0, Step 29 is bypassed. 

til 

(29) If J is not equal to zero, the J instrument state is set 
equal to 1. 

(30) If S equals S', no change in failure status has occurred, 
and the correction matrix calculation is bypassed. 

(31) The new B matrix is calculated from Eq. 3-14 with the 
rows of A corresponding to the failed instruments taken 
as zero. See Steps 11 to 17 in Paragraph 3.2.2 for the 
detailed description. 

(32) The package input is estimated from the unfiltered instru- 
ment outputs. 

3.2.4 Bayesian Decision Theory Algorithm 

The Bayesian Decision Theory Algorithm is described in Refer- 
ences 15 and 16. Actually, there are two algorithms derived in 
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Reference 15: one for "mean failures" and one for "variance failures." 
Apparently these terms refer to the extremes of constant and white noise 
errors. Reference 15 is not completely transparent to this reader, and 
since Reference 16 sets forth the mean failure algorithm clearly, the 
mean failure algorithm will be chosen. 

The mean failure algorithm presented in References 15 and 16 is 
in error, in that it produces incorrect decisions for failures that lead to 
errors with negative polarity. This shortcoming can be remedied easily 
by the judicious introduction of absolute value operations into all of the 
decision equations. 

Once this is done, it becomes obvious that the 21 decision equations 
set forth for six instruments involve only six different test signals and 
that these test signals are proportional to the residuals, e, of Eq. 3-17. 
Furthermore, the decision tree reduces to the selection of the test sig- 
nal having the largest absolute value and the comparison of it to a thresh- 
old. If it exceeds the threshold, the associated instrument is failed. 

If it does not, zero instruments are failed. 

When one instrument is failed, the decision equations for the remain- 
ing five instruments must be implemented. The test signals used in 
References 15 and 16 are arbitrary and have unequal variances represent- 
ing the errors of the other instruments^’in unfailed operation. Considera- 
tion of the symmetry of the dodecahedron configuration leads one to 
believe that, if the test signals were chosen to minimize their variances, 
their variances would all be equal to each other and less than or equal to 
the variances of the test signals of References 15 and 16. 

This guess can be shown to be correct, and the minimum variance 
test signals turn out to be the residuals, 7 , of Eq. 3-23. Therefore they 
will be used in place of the original test signals. Again the decision tree 
reduces to the selection of the test signal having the largest absolute 
value, and the comparison of it to a threshold, to determine whether or 
not a second failure is present. 

To reduce the frequency of false alarms, Reference 16 calls for two 
consecutive determinations of a failure before an instrument is switched 
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out. However, this process will cause an increase in missed alarms 
not taken into account in the previous analysis in Reference 16. It seems 
clear that this process is approximately equivalent to raising the detection 
threshold. Therefore, it will be ignored. The thresholds will be set 
experimentally rather than theoretically at any rate. 

Reference 16 also calls for two sampling periods, 10 min for soft 
failures and 0.01 sec for hard failures. The 10-m.in period samples are 
averaged over the preceding sampling period. This procedure represents 
another prefilter, resembling somewhat the SIRU prefilter. The simple 
low-pass filter will be used instead for the same reasons as given for the 
Squared Error Algorithm in Paragraph 3. 2. 3. 

It is clear from the foregoing that considerable liberties have been 
taken with Gully's original concept. This course of action was necessary 
because his scheme was further from a "reduction to practice" than the 
others. In summary, when no failures have been detected, the largest 
residual of Eq. 3-17 is compared against a threshold. If it exceeds the 
threshold, the corresponding instrument is failed. When one failure has 
been detected the largest residual of Eq. 3-23 is compared against a 
threshold. If it exceeds the threshold, the corresponding instrument is 
failed. No healing is permitted. Figure 3-13 shows the algorithm flow 
diagram. 

The Bayesian Decision Theory Algorithm has three adjustable param- 
eters: the prefilter time constant and the threshold levels TH^ and TH^. 

The flow diagram description follows: 

(1) The filter 3 vector f, the instrument state S, the failed 
status indicator k, and the correction matrix B are initi- 
alized before the trial. 

(2) Each minor cycle the failed status indicator is tested. If 
it equals 7, two failures have occurred, and the failure 
detection and isolation logic is bypassed. 

(3) If 0 or 1 failure has occurred, the inconsistency state 
vector v is calculated and filtered. 

(4) If k is not equal to 0, one failure has been isolated, and 
the first failure detection logic is bypassed. 
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(5) If no failures have been isolated, the residuals are 
calculated, and the loop to find the maximum residual is 
initialized, j is the number of the instrument having the 
largest residual and u is the value of its magnitude. 

(6) If the magnitude of the i^ 1 residual is less than u max > then 
Step 7 is bypassed. 

(7) j is set equal to i, and Uj^x is set equal to the magnitude 
of u. if | u. | is greater than u max - 

(8) If the magnitude of the largest residual is less than the 
threshold, no failure is detected and the rest of the fail- 
ure detection logic is bypassed. 

th 

(9) Instrument j is failed, k is set equal to j. The j instru- 
ment state is set equal to 1. m is set equal to 1 to show 
that a change has occurred in the instrument states S. 

The matrix for calculating the residuals from the filter 
vector when instrument k is failed is calculated. 

(10) If one failure has already been isolated, m is set equal to 
0 to show that no change in S has occurred. 

(11) The residuals (excluding instrument k) are calculated. The 
loop to find the largest residual is initialized. 

(12), (13) The largest residual is found as in Steps 6 and 7. 

(14) If the magnitude of the largest residual is less than the 
threshold, no second failure is detected. 

(15) If a second failure is detected, k is set equal to 7 and the 
jth instrument state is set equal to 1. 

(16) If no second failure is detected, the correction matrix 
calculations are bypassed unless m equals 1, showing 
that a first failure was detected in this minor cycle. 

(17) The new B matrix is calculated from Eq. 3-14 with the 
rows of A corresponding to the failed instruments taken 
as zero. See Steps 11 through 17 in Paragraph 3.2.2 for 
the detailed description. 

(18) The package input is estimated from the unfiltered instru- 
ment outputs. 
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Figure 3-13. Bayesian Decision 
Theory Algorithm 
Flow Diagram 
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3.2.5 Maximum Likelihood Algorithm. 


The Maximum Likelihood Algorithm is described in References 17 
and 18. A maximum likelihood technique is used to obtain estimates of 
the instrument variances from the inconsistency vector v. The instru- 
ment variances are then used in a weighted least-squares estimator. 

The inconsistency vector is modified by subtraction of the residuals to 
reduce the interaction of successive failures. The squared residuals 
are added to the variance estimates to compensate for the decrease that 
would otherwise result from the subtraction of the residuals. 

The prefilter in the original algorithm, which was imbedded within 
the algorithm itself, proved to be unsuccessful. Therefore the instru- 
ment outputs are prefiltered prior to use in the algorithm. Figure 3-14 
shows the algorithm flow diagram. The Maximum Likelihood Algorithm 
has two adjustable parameters, the prefilter time constant and the a 
priori error variance Q q . 

The flow diagram description follows: 

(1) The residuals e., the _filter 6 vector f, the previous value of 
the weighting mktrix Q, and the correction matrix B are 
initialized before the trial. 

(2) Each minor cycle the filter outputs are calculated. The 
residuals are subtracted from the filter outputs, and the 
inconsistency state vector is calculated. The maximum 
likelihood variance estimates Q' are found and modified 
to account for a priori information and to compensate for 
the subtraction of the residuals to give the instrument 
covariance matrix Q. The iteration counter is initialized. 

(3) If the change in any one of the diagonal elements of the 
covariance matrix is greater than the factor CRIT times 
its previous value, a new correction matrix is required. 

(4) The iteration count is increased by one. The previous 
value of the covariance matrix is set equal to the new 
value. The correction matrix is updated. 

(5) The new residuals are found. The new covariance matrix 
is obtained, (v is now zero, so that Q' is eliminated from 
the expression for Q. ) 

(6) If 100 iterations have occurred, the iteration loop is 
terminated. 
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(7) If the change in any one of the diagonal elements of the 
covariance matrix is greater than the factor CRIT times 
its previous value, a new correction matrix is required. 

(8) The package input is estimated from the unfiltered 
instrument outputs. 

3. 2. 6 Minimax Algorithm 

The Minimax Algorithm is described in References 19, 20, 21, 22, 

23, 24, 25, and 26. It uses the 15 parity signals described in Para- 
graph 3. 1. 2. 3. The original algorithm uses no filtering at all, so a pre- 
filter in the form of a simple low-pass filter has been provided to put it 
on an equal footing with the other algorithms. This prefilter rejects the 
large, high-frequency rate errors caused by the instrument quantization. 

The original algorithm (as detailed on a MAC listing) solved for the 
parity equation weights (±c and ±s in Eq. 3-5) each time through in order 
to handle the gimballed case where the relative instrument orientations 
can change. As pointed out by the authors, this procedure is not necessary 
for the dodecahedron strapdown case. Thus the formulation of Eqs. 3-57 
and 3-58 is used instead. The parity signals used in the Minimax Algo- 
rithm are normalized by a factor of l/(2c+2s) so that, if the instrument 
outputs take on the values ±1, the magnitude of the normalized parity 

signal lies in the range ±1. In FALLSIM this factor is included in the con- 

2 

version of the threshold levels from pg and deg/hr to m/sec and rad/sec. 

Failure detection takes place when any one of the parity signals 
exceeds a threshold. The original algorithm has a technique for finding 
the maximum parity signal without calculating all of them. When detec- 
tion occurs, a second loop calculates all of the parity signals to perform 
isolation. Since the calculation of the parity originals is much simplified 
by the use of Eqs. 3-57 and 3-58, this procedure is modified. Both 
detection and isolation are combined in a single loop in which all of the 
parity signals are calculated. Neither this change nor the one in the pre- 
ceding paragraph changes the action of the algorithm in any way. 

For isolation of the first failure, it is assumed that only one failure 
can occur at a time. The four instruments that contribute to the value of 
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a particular parity signal are called a quartet. If the signal exceeds its 
threshold, the quartet is dirty; if not, it is clean. If no instruments 
have failed, all quartets must be clean. Thus, a failure is detected if 
any quartet is found to be dirty. 

With six instruments and one failure, there must be one set of five 
good instruments. With a set of five instruments called a quintet, the 
isolation algorithm finds every quintet whose quartets are all clean. 

Such a quintet may be a set of good instruments and will be called a clean 
quintet. An instrument which is excluded from every clean quintet must 
be failed since a good instrument would be in some clean quintet. Thus 
the failed instrument is the one excluded from all clean quintets. 

For the second failure, there are five parity signals which are inde- 
pendent of the output of the first failed instrument. When four of these 
exceed their thresholds, the second failed instrument is successfully iso- 
lated. The zeros in the equation for the sole clean parity signal corre- 
spond to the two failed instruments. The algorithm can isolate those 
double simultaneous failures which cause all but one parity signal to 
exceed its threshold. 

It turns out that if a first failure causes all of the parity signals to 
which it contributes with a coefficient of ±c to exceed the threshold, the 
Minimax logic will cause the failure to be isolated. However, in order 
to be isolated, a second failure must cause four parity signals to exceed 
the threshold. It contributes to two of them with a coefficient of ±c and 
two with a coefficient of ±s. Thus, in effect, the threshold for second 
failures is higher than the threshold for first failures by the ratio c/s 
or 4. 2 dB. 

Failure correction may be accomplished either with the least- squares 
technique or with the Bounding Sphere Algorithm. This algorithm mini- 
mizes the maximum possible estimation error. It is capable of modify- 
ing its action when failures are detected but not isolated unlike the least- 
squares algorithm. See the references for a description of this algorithm. 

Figure 3-15 shows the algorithm flow diagram. (The Bounding 
Sphere Algorithm flow diagram is not included since none has been seen 
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by the author. The Bounding Sphere FORTRAN coding was obtained by 
translating MAC. ) The Minimax Algorithm has two adjustable parameters, 
the prefilter time constant t^, and the threshold level THM. The Bounding 
Sphere Algorithm also requires a threshold level THS. In the absence 
of a prefilter, these levels are theoretically identical. The original 
algorithm makes no distinction between them. 

The flow diagram description follows: 

(1) The filter 3 vector, the list of the numbers of the active 
instruments MNj, the number of active instruments n, 
the number of detected but unisolated instruments CNO, 
the instrument state S, and the correction matrix B are 
initialized before the trial. 

(2) Each minor cycle the inconsistency state vector v is cal- 
culated and filtered. The present value of S is stored. 
QUINTj, representing the quintet of instruments excluding 
instrument i, is zeroed. IFAIL and NUM are set equal 

to 0. The parity signal index p is set equal to 1. 

The four indices i, j,k, and 1 are the numbers of the instruments 
included in the p quartet. The four DO loops cause p to run from 1 to 
15 selecting the rows of the matrix F in Eq. 3-58. Thus the parity sig- 
nals of Eq. 3-53 are calculated in the order shown in Eq. 3-54. 

th. 

(3) If none of the instruments in the p quartet are failed, IT 
is 0. 

(4) If the quartet includes a failed instrument, it is ignored. 

tVi 

(5) The p parity signal is calculated. 

(6) The parity signal is tested against the threshold. 

(7) If the signal is less than or equal to the threshold, the 

p fc k quartet is clean and the number of good quartets, NUM, 
is increased by 1. 

(8) If NUM is not equal to 1, Step 9 is bypassed. 

(9) The numbers of the instruments in the first good quartet 
are stored. (They will provide the solution for a double 
simultaneous failure if there is only one good quartet, 

i. e. , if NUM = 1 in Step 29. ) 
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(10) - (13) If the signal is greater than the threshold, the 

pth quartet is dirty. The m** 1 quintet excludes instrument 
m. Therefore if i, j, k, and 1 are all not equal to m, the 
m^ 1 quintet includes the p“* quartet and must be dirty, 
since its quartets are not all clean. If i, j, k, or 1 equals 
m, the quintet does not include one of the instruments of 
the quartet, so that the p th quartet says nothing about the 
quintet. 

(14) If the m^k quintet is dirty, QUINT^ is set equal to 1. 

(15) IFAIL is set equal to 1 to show that at least one parity 
signal has exceeded the threshold. 

(16) The parity signal index is increased by 1. 

(17) If IFAIL equals 0, no failures have been detected. 

(18) CNO is set equal to 0. 

(19) If n equals 4, only 4 instruments are unfailed. The new 
failure indicated by IFAIL nonzero cannot be isolated. 

(20) The number of detected but unisolated failures is assumed 
to be 1, CNO is set equal to 1. 

(21) m is set equal to -1. 

(22) If the i^ quintet is dirty, QUINT^ is nonzero, and the rest 
of the loop is bypassed. 

(23) If m is not equal to -1, this is not the first clean quintet. 

(24) CS is set equal to the number of the first clean quintet. 

(If it turns out to be the only clean quintet, then CS will 
be the number of the failed instrument. ) 

(25) The clean quintet counter is increased by 1. 

(26) If m is less than 0, there were no clean quintets. If m 
equals 0, there was one clean quintet. If m is greater 
than 0, there were two or more clean quintets. 

(27) If NUM = 0, there were no clean quartets. 

(28) If there were no clean quartets, it is assumed that there 
are three good instruments remaining. The number of 
failures detected but not isolated is the number of active 
instruments less 3. 

(29) If NUM = 1, there was one clean quartet. 
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(30) The number of active instruments, m, is set equal to 4. 

The list of the numbers of the active instruments is 
changed to the values stored in Step 9. The number of 
detected but unisolated failures is 0. 

(31) If there were two or more clean quartets, then 2 detected 
but isolated failures are assumed. 

(32) There was one clean quintet. The index k is set equal 
to 1. 

(33) If MN. equals CS it is the number of the failed instrument. 

(34) The new list of numbers of the active instruments is found 
excluding the failed instrument, k is increased by 1. 

(35) The number of active instruments is reduced by 1. The 
number of detected but unisolated instruments is set 
equal to 0. 

(36) The instrument state is updated. 

(37) There were two or more clean quintets. One detected but 
unisolated failure is assumed. 

(38) If LSQ equals 0, the Bounding Sphere Algorithm is to be 
used. If LSQ is not equal to 0, the least-squares algo- 
rithm is to be used. 

(39) If S equals S', no change in failure status has occurred, 
and the correction matrix calculation is bypassed. 

(40) The new B matrix is calculated from Eq. 3-14 with the 
rows of A corresponding to the failed instruments taken 
as zero. See Steps 11 through 17 in Paragraph 3.2.2 
for the detailed description. 

(41) The package input is estimated from the unfiltered instru- 
ment outputs. 

(42) The Bounding Sphere Algorithm estimates the package 
input making use of CNO, which is ignored by the least- 
squares algorithm. 

3. 2. 7 Adaptive 72 Algorithm 

The Adaptive 72 Algorithm is described in References 27 and 28. 

It is the most complex of the eight algorithms studied and one of the most 
recent. Unfortunately, these two features combined to make it impossible 
to include it in FAILSIM as planned. The complexity mandated the use of 
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the detailed program flow diagrams in writing the routines for FAILSIM 
but the diagrams were unavailable in time. 

V 8 Sftrmential Algorithm 

’ ' — _ - ~ ~i ,.r 7—; a— — 

The Sequential Algorithm is described in Reference 29. It is also 
complex and new. The implementation presented here is based upon the 
algorithm as it existed when the reference was written. It has not been 
possible to include in FAILSIM the subsequent improvements made by its 
authors. 

The Sequential Algorithm consists of a modified Kalman-Bucy filter. 
A detailed description of the filter is given in Paragraph 3. 1.2.4. In 
this section we shall see how it is modified to perform FDIC. 

In the original algorithm, any instrument whose error state exceeded 
a threshold was taken as failed. This logic has been modified so that only 
the instrument whose error state exceeds the threshold by the greatest 
amount in any single pass through the filter equations is taken as failed. 

The nonlinearity required for failure detection is provided by two 
distinct threshold levels. One level is set quite high and is applied to 
the Kalman-Bucy residuals. If the threshold is exceeded, the associated 
instrument is excluded from the measurements. When the residual is 
again less than the threshold the instrument is reinstated. The second, 
much lower level, is applied to the instrument error states. An instru- 
ment whose error state exceeds the threshold by more than any other 
error state during a single pass through the filter equations is excluded 
permanently from use. The residual threshold is constant. The error 
state threshold is proportional to the standard deviation of the error state, 
as given by the filter covariance matrix. Figure 3-16 shows the algorithm 
flow diagram. The Sequential Algorithm has nine adjustable parameters 
as follows: 

(1) Initial variance for package input states 

>!< 

(2) Initial variance for instrument error states 

(3) State noise variance for package input states 


The standard deviation is input to the program. 
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Initialization 



FQTDOUT frame 


Figure 3-16. Sequential Algorithm 
Flow Diagram 
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(4) State noise variance for instrument error states 

* 

(5) Variance for measurement errors 

(6) Time constant for package input states 

(7) Time constant for instrument error states 

(8) Threshold for residuals 

(9) Ratio of threshold for instrument error states to instru- 
ment error standard deviation. 

The flow diagram description follows: 

(1) The estimated state 9 vector x, the state covariance matrix 
P, the state transition matrix $ , the state noise covari- 
ance matrix Q, the measurement noise matrix R, and 

the measurement matrix H are initialized before the trial. 

(2) Each minor cycle the state vector and covariance matrix 
are advanced. 

(3) Temporary working values of the state vector and covari- 
ance matrix, x' and P’, are stored. The instrument state 
is zeroed. 

(4) The j** 1 residual is calculated. 

(5) The j* residual is tested against the high level threshold. 

(6) If the residual is less than the threshold, the filter matrix 

K. is calculated, and x 1 and P' are reset. 

J 

(7) If the residual is greater than the threshold, the j instru- 
ment state is set equal to 1. 

(8) The low-level failure detection and isolation process 
begins by setting i and SMAX to 0. 

(9) If the j^ 1 instrument was excluded from the current reset, 
IT will be nonzero. 

th 

(10) If the j instrument was excluded from the current reset, 
low-level failure detection is not performed. 

th 

(11) The estimate of the j instrument error is squared. 


*The standard deviation is input to the program. 
**Used to calculate the state transition matrix. 
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(12) If the square is less than the squared threshold times the 
instrument error variance, the instrument is not failed. 

(13) - (14) At the end of the loop, SMAX will equal the squared 
error of the instrument having the largest error, and i 
will equal the number of the same instrument. If the error 
of no instrument exceeds the threshold, SMAX and i will 
be equal to 0. 

(15) If i is not equal to 0, a failure has been detected. 

th 

(16) If a failure has been detected, the i instrument error 
estimate is set to a large value. Subsequently, the state 
of that instrument will always be set equal to 1 in Step 7. 
The filter calculations are now repeated with instrument 
i excluded. 

(17) If no failures have been detected, the state vector and 
covariance matrix are updated from the working values. 

The first three elements of x comprise the output to the 
strapdown algorithm. 
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4. STRAPDOWN INERTIAL PACKAGE MODEL 


The strapdown inertial package model simulates the operation of 
the package in the navigate mode only. The calibration and alignment 
functions are assumed to have been completed prior to the commencement 
of a computer run. Thus the errors simulated do not represent the total 
instrument and package errors, but only the residual errors remaining 
after calibration. 


Because of computing time restraints, it is out of the question to 
simulate a cruise-type system for any significant length of time. Thus a 
boost-type system has been chosen. It should be adequate to allow the 
candidate algorithms to compete meaningfully against each other. 

4. 1 INSTRUMENT CONFIGURATION 


The gyro and accelerometer input axis orientations with respect to 
the package axes are defined by the A matrix given in Eq. 3-2. It is 
assumed that the gyros and accelerometers are defined by the same A 
matrix (that is by the same dodecahedron orientation) because, although 
it is conceivable that a system could be designed with different dodecahe- 
dron orientations for the gyros and the accelerometers, nobody has done 
so, so that the added complexity required in the simulation to treat this 
case would probably be wasteful. 


The orientation of each gyro and accelerometer about its input axis 
is defined by the angles 0g^ - 0g^ and 0 a ^ - 0 a ^> are positive 

rotations about the input axis from a reference orientation. When an 
instrument is in the reference orientation, its output axis lies along the 
positive direction of one of the package axes. From the matrix A and the 
12 orientation angles, it is possible to determine the 12 matrices 
AG1-AG6 and AA1-AA6 that relate individual gyro input, output, and spin 
axes and accelerometer input, pendulous, and output axes to package axes. 
We have, for the gyros. 


X input 

X 

output 

X spin 


= AGi 


'Package 


(4-1) 
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and, for the accelerometers 
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where 


(4-9) 


c. 

= cos 0 or 

cos 0 

1 

*i 

a. 

l 

s. 

= sin 0 or 

sin 0 

1 

*i 

a i 


til 

Note that the first row of AGi and AAi equals the i row of A. 

Except for a relabeling of instruments and a rotation of 180 deg 
about certain of the instrument axes, the reference instrument orientations 
correspond to those of Gilmore 1 s SIRU. The instrument labels and 
180 deg rotations have no significance to a simulation of this type. 

Table 4-1 shows the relabeling required and the axes about which 
a 180-deg rotation is required to achieve correspondence between the! 
gyro axes for the reference orientation and SIRU. 

In SIRU the accelerometer input and output axes are the same as for 
the gyros, and the pendulous axis is the negative of the gyro spin axis. 

The same relationship holds in the simulation whenever the 0g and 0 a 
are identical. 1 

The MSFC breadboard dodecahedron (BB DDH) uses single-axis 
platforms, so that the directions of the output and spin axes vary with 
time. This feature is not modeled in the simulation, as it is felt that it 
would have no impact on the relative performance of the FDIC algorithms 
and would increase computer time requirements. 


Table 4-L Simulation — SIRU Correspondence 


Simulation 

SIRU 

Gyro 

Gyro 

Axis 

1 

C 

Output 

2 

D 

Spin 

3 

E 

None 

4 

F 

Spin 

5 

N A 

None ' 

6 

B 

Input ' 


This rotation may be 
accomplished by 0g^ = 


180 deg. 
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The package axes are related to the body axes by a matrix R 


x , = R x, , 

package body 


(4-11) 


where R is defined as the product of three single-axis orthogonal rotation 
matrices 


R = T 


(V S) T (V 1r 2 ) t (V iR i) 


th 


(4-12) 


The 


where T(x, m) is a rotation through the angle x about the m axis. 

angles Or , Or , and Or and the axis indicators 1r , 1r , and i R are 
1 2 3 1 2 3 

all read in as data. Note that the rotations from body-to-package axes 

are about axes iR^> iR^> an< ^ ^R^ through angles 9^, ant ^ ^R^’ 

respectively, in that order. Since iR^, iR > and 1 r^ can each take on 

any of the values 1, 2, or 3, any order of rotations is possible. 

4. 2 UNFAILED INSTRUMENT ERRORS 


The unfailed instrument simulation must exhibit errors, since it is 
the errors of the normally operating instruments that make the failure 
detection process difficult. However, the simulation of errors should not 
be excessively detailed and exhaustive because of the resultant penalties 
in computer time and storage. Therefore, a limited set of relatively 
easily calculated errors has been selected. Their combined character- 
istics are complex enough so that they do not correspond to any of the 
simple error models hypothesized in the derivation of the different 
algorithms. For the numerical values of the errors, see Paragraph 6. 2. 4, 
Appendix A, and Reference 30. 


4. 2. 1 Gyro Errors 

The error in a gyro output is given by 



where the errors are: 


(1) Bias drift rate 

(2) Random drift rate 

(3) Scale factor error 

(4) Input axis misalignment about spin axis towards output 
axis 

(5) Input axis misalignment about output axis towards spin 
axis 

(6) Acceleration sensitivity along input axis (spin axis 
mass unbalance) 

(7) Acceleration sensitivity along output axis 

(8) Acceleration sensitivity along spin axis (input axis 
mass unbalance) 

(9) Anisoelastic drift rate. 

The standard deviations of the quantities E^ through E^ are input. 
Independent values of E^ through E^ are obtained for each gyro by the 
Monte Carlo method. 

In addition to the above errors, there are errors introduced by the 
sampling and quantization process. The gyro output angular velocities 
are integrated over one sampling period. The result is then quantized. 
The error resulting from the quantization process is retained and added 
to the result of the integration over the next sampling period. This tech- 
nique simulates the storage of quantization error information in the gyro 
float that occurs in an actual gyro. 

4. 2. 2 Accelerometer Errors 

The error in an accelerometer output is given by 

c = E. + E_ n + E_ a, + E . a + E_ a 
1 2 3i 4 p 5o 

where the errors are 

(1) Bias 

(2) Random error 
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(3) Scale factor error 

(4) Input axis misalignment about output axis towards 
pendulous axis (or pendulous axis acceleration 
sensitivity) 

(5) Input axis misalignment about pendulous axis towards 
output axis (or output axis acceleration sensitivity). 

The quantities E^ through Ej. are input as for gyros. The acceler- 
ometer outputs are integrated, sampled, and quantized in the same manner 
as the gyro outputs. 

4. 3 STRAPDOWN ALGORITHM 

The incremental velocity and angle outputs of each FDIC algorithm 
go to a strapdown algorithm that performs the attitude and inertial velo- 
city calculations. All of the strapdown algorithms are identical. 

The gravity calculations have only a minor effect on the propagation 
of errors in a boost inertial navigation system. Therefore, they are 
eliminated from the trajectory generator and each strapdown algorithm 
in order to save computer time and storage. The velocity errors at injec- 
tion have a much stronger effect on navigational performance than do the 
position errors. Since position is no longer needed for the gravity calcu- 
lations, the position integration can also be eliminated from the trajectory 
generator and each strapdown algorithm. These simplifications should 
have very little effect on the performance of the FDIC algorithms. 

The strapdown algorithm is based on the approach described in 

Reference 31, with an improved attitude algorithm. The attitude algorithm 

32 

chosen is McKern's third-order quaternion algorithm that is used in 
SIRU. This algorithm is accurate and efficient and includes a correction 
for computation errors. The incremental quaternion Ap is calculated 
from the incremental angle vector <j> by 

- 1 - HMi - 4 W ± )* - <4 - 14) 

where 4> * is the incremental angle vdctor from the' preceding sampling 
period. The quaternion representing the vehicle orientation is updated by 

p — p Ap (4-15) 
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33 34 

The velocity in body axes is updated by Yachter's algorithm ‘ 

V B~ V B + Av - <I>X ( V B + Tt) (4_16) 

At a slower frequency, Vg is transformed into inertial axes and the result 
is added to the inertial velocity v^. At the same time, Vg is zeroed in 
Eq. 4-16. 

Vj— Vj + B Vg (4-17) 


where B is defined in terms of p 

2(p 3 p 1 + p 2 p 4 ) 

p/ 2(p 2 P 3 - p l p 4> 

2 2 2 2 

- Pl ~ p 2 + p 3 + P 4_ 

(4-18) 

4.4 TRAJECTORY 

The trajectory consists of two parts, a nominal trajectory and 
superimposed vibrational motions. The trajectory remains the same 
throughout a run. The statistics of the vibrational motions remain the 
same throughout a case. The actual vibrational motions are Monte Carloed 
and vary from trial to trial. 

4. 4. 1 Nominal Trajectory 

The nominal trajectory is specified to FAILSIM by means of a sub- 
routine, TRAJ. The acceleration and angular velocity are calculated in 
TRAJ by piecewise analytical functions of time which are written in 
FORTRAN. Discontinuities are permitted between the pieces, provided 
they occur at integral multiples of the minor cycle period At. 

4. 4. 2 Vibrational Motions " 

The translational and rotational vibrational accelerations are 
generated by means of gaussian pseudorandom numbers, which are filtered 


p l 2 - p 2 2 - p 3 2 + p 4 2 2(p l p 2 - p 3 p 4> 


2 ( p l P 2 + P 3 P 4 ) 


2( p 3 P x - P 2 P 4 ) 


2 2 2 , 
-Pl + Po - P* + 


2 ( p 2 P 3 + P x P 4 ) 
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by first-order difference equations so as to obtain first-order Gauss- 
Markov random processes. If unmodified, such accelerations would 
cause the translational and rotational velocities to random walk away 
from the nominal. To prevent this, the accelerations are passed through 
high-pass filters made by feeding back the velocity and position changes 
caused by the vibration through gains as shown in Figure 4-1. In the 
figure, n is white noise, m is the Gauss-Markov process, and a, v, 
and p are the translational and rotational acceleration, velocity, and 
position, respectively. In the cases of the three rotational vibrations 
and two of the translational vibrations, the feedback can be considered 
to represent the vehicle attitude control and guidance systems. For 
the vehicle longitudinal axis, no physical significance can be ascribed 
to the translational vibration feedback, and the gains may be set to 
zero if the random walk can be tolerated in this axis. 



Figure 4-1. Vibrational Motion Generator 

The first-order difference equation which represents the low-pass 
filter that generates the Gauss-Markov process is \ 

m. , . = K m. + K, n. . (4-19) 

l+l a i b i+ 1 
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In order to make m correlated with time constant t, we require that 


<\ m i+l m i> 

(-i 2 ) “ 

where T is the sampling period. Therefore, 

-T 

K = e T 
a 

The constant will be adjusted to give the desired variance of the 
vibrational acceleration at the output of the high-pass filter. First let us 
find the Z transform of Eq. 4-19 

z m = K m + K, z n (4-22) 

a b 


(4-20) 


(4-21) 


Then 


G.(z) = — 
1 ' m 



(4-23) 


Now let us find a sampled data approximation to the high-pass filter from 
Figure 4- 1 , 


°2 


2 

s 


s 2 + K s + K 
v p 


(4-24) 


Let 


~ 2 z - 1 
“ T z + 1 


(4-25) 


then 


2 

B z + B . z + B ? 

G 2 ( z ) = — ~2 " 1 ( 4 - 26 ) 

z + A^ z + 
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where 


B 

o 



/ 2 \ 

2 \ K P T - 4 i 


4+2 

K 

T + K 

T 2 


V 

P 


4 - 2 

K 

T + K 

T 2 


V 

P 


4 + 2 

K 

T + K 

T 2 


V 

P 




4 


4 + 2 

K 

T + K 

T 2 


V 

P 



(4-27) 


(4-28) 

(4-29) 


B. = -2 B 
1 o 


(4-30) 


The combined transfer function is 


K,z B z + B. z + B- 

G(z)=G 1 (z)G 2 (z)= r V* ~ °2 ~~~ 7~ ~T " ~ 

a z + Aj z + A 2 


(4-31) 


The power spectral density of the random number sequence is 


$ (z) = 

nn ' ' 


2ttT 


(4-32) 


and the power spectral density of the filter output is 


S (z) = G (z) G (z" 1 ) ® (z) 

aa nn 


(4-33) 


The variance of the filter output is 


a 2 = — (f) $ (z) z - ^ dz (unit circle) 

a j Jr aa 


(4-34) 
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From Eqs. 4-31, 4-32, 4»33, and 4-28 


2 K b 2 I 1 

i B z 
( o 

2 + B x z + B 2 ) ( 

' _2 _ i \ 

B z + B.z + B- 
i o 1 Z ) 

z _1 

5 ^ J (- K a)l 

/ 2 ' 

z 4 + A. z + A, 

^ 1 2 , 

) ( z“ 2 + A i z_1 + A 2 

1' 

( z_1 - K a 


(4-35) 


From the residue theorem 



2 Res 



(4-36) 


where Pj^ are those poles of the integral which lie within the unit circle; 

thus K, is 
b 


K b = 


V rRes ( P k) 


(4-37) 


The residues are evaluated, and is calculated in subroutine FKB. 

For the angular vibrations, the velocity as well as the acceleration 
is required. The transfer function from m to v is 


G 3 (s) 


s 2 + K s + K 
v p 


If we apply Eq. 4-25 again. 


G 3 < Z) 


C z 2 + C.z + C, 
o 12 

2 

z + A^z + A 2 


(4-38) 


(4-39) 


where 


C 

o 



2T 

j 

4 + 2K T + K T 
P P 


= 0 


= -C 

o 


(4-40) 

(4-41) 

(4-42) 
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Figure 4-2 shows a block diagram of a sampled data system which 
can give both a and v. The equations which must be solved by the computer 
are 


m *- K m 4- K, n 
a b 

(4-43) 

a l — a 2 

(4-44) 

a 2 a 3 

(4-45) 

a 3 ■*" m - A x a 2 - A 2 &1 

(4-46) 

a *" B o a 3 + B 1 a 2 + B 2 a l 

(4-47) 


v *’ — C a 0 + a. 


(4-48) 















in that order. Initially a^, a 2> and a 3 are set equal to zero; m is initially 

K,n 

m = b — (4-49) 

The translational acceleration experienced by the strapdown package 
is a combination of the translational acceleration of the center of mass of 
the vehicle and the effect of the rotational acceleration and velocity on the 
lever arm between the center of mass and the package. This effect may 
be represented by 


a TOT 1 = a TRANS 1 ~ R | V ROT 2 + v ROT^j 

(4-50) 

a TOT 2 = a TRANS 2 + R a ROT 3 

(4-51) 

a TOT 3 = a TRANS 3 “ R a ROT 2 

(4-52) 


where the lever arm, of length R, is assumed to lie along the 1 axis. 

4. 5 USE OF INPUTS FROM REAL INSTRUMENTS 

The simulation program is capable of accepting the outputs of real 
instruments on a strapdown package mounted on a test table and bypassing 
the computation of simulated instrument outputs. Since the MSFC package 
has only three accelerometers, an option is provided to convert their 
three outputs into six outputs, which the FDIC then converts back into 
three outputs. Since the simulation program has no calibration, alignment, 
or compensation calculations, the outputs of real instruments must be suit- 
ably compensated before being input to the program. Failures are to be 
induced in the instrument hardware. 

A strapdown package on a test table is unaccelerated, so that the 
ideal inertial velocity can be calculated without knowledge of the package 
orientation. Since the FAILS IM strapdown algorithms have no gravity 
calculation, the ideal inertial velocity cannot be obtained merely by trans- 
forming the test table velocity into inertial axes. Instead, the total sensed 
acceleration must be integrated. 
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The initial inertial velocity may be taken as zero, and the initial 
acceleration of gravity vector may be taken as lying in the X-Z plane of 
the inertial coordinate system without loss of generality. Thus the acceler- 


I cos L cos t | 

cos L sin t > (4-53) 

sin L f 

and the velocity vector is the integral of Eq. 4-53 

^ cos L sin £2 t 

^ cos L. (1 - cosfit) 
t sin L 

where 

L = astronomic latitude of test table 
g = magnitude of gravity vector at test table 
£2 = earth rate. 

The initial package orientation must be known in order to initialize 
the quaternion in each strapdown algorithm. The quaternion defining 
up-east-north (UEN) axes with respect to inertial axes is 

= cos -^r - j sin (4-55) 




The initial orientation of the package axes with respect to UEN axes is 
represented by rotations about three axes taken in any order. The total 
quaternion describing the initial orientation of the body axes with respect 
to inertial axes is 

P = P L P^P^ ipj) P(&P 2 > iPz) P(0P 3 . iP 3 ) (4-56) 
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where p(x,m) is a quaternion representing a rotation through the angle x 
th 

about the m l axis. The angles 0p^, 9p£> an ^ ^Pj anc ^ the axis indicators 
ipj, ip^, and ip^ are all read in as data. Note that the rotations from 
UEN to package axes are about package axes ip^, ip 2 » and ip^, through 
angles 0p^, and Sp^, respectively, in that order. Since ip^, ip 2 » 

and ip^ can each take on any of the values 1, 2, or 3, any order of 
rotation is possible. 

Since FAILSIM has no way of knowing the test table rotations, the 
ideal quaternion computation cannot be performed. Therefore, the test 
table should be returned to its initial orientation at the end of the run so 
that the initial quaternion can also serve as the ideal quaternion for cal- 
culation of the orientation errors at the end of the run. 
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5. FAILURE MODES 


A representative set of gyro and accelerometer failure modes is 
chosen. Although all possible failure modes are not included, and those 
included may be represented in a somewhat simplified manner, they are 
sufficiently varied in nature to give the competing FDIC algorithms a 
good testing. Both hard and soft errors are included. 

There are two types of failure modes to be simulated. The additive 
failure modes add to the output produced by the instrument in unfailed 
operation. For example, a shift in the bias drift would not alter the 
response of a gyro to inputs. The substitutional failure modes substitute 
an incorrect output for the correct one. The incorrect output is indepen- 
dent of the output that would be produced by the instrument if it were 
unfailed. For example, a failed instrument could produce zero output. 

Each failure will be characterized by the following parameters: 

• Time of failure 

o Time failure disappears (for representing glitches) 

« Which instrument fails 
9 Failure mode. 

In addition, certain failures will be characterized by algebraic values 
or standard deviations, as described below. 

5.1 GYRO FAILURES 

5. 1. 1 Additive Failure Modes 

Mode 1: Bias Drift Rate Shift 

The gyro drift rate shift is a step function in angular velocity. 

Either the algebraic value of the shift may be input directly or it may be 
Monte Carloed. If it is Monte Carloed, the amount of the shift is obtained 
as a gaussian pseudorandom number with a specified standard deviation. 

The gaussian model for the underlying process that causes shifts 
in bias drift is not based on any actual failure data, but seems to be a 
reasonable choice. However, for there to be a failure, the shift must be 
greater than some threshold value. This value must be specified before 
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one can accumulate empirical failure rate or MTBF data on a real gyro. 
Obviously, the lower this value is, the higher the failure rate will seem 
to be. This threshold value is input to the program along with the stan- 
dard deviation. If the magnitude of the gaussian pseudorandom number is 
less than the threshold value, it is rejected and a new random number is 
chosen. 

The threshold value should be chosen high enough to exclude '’fail- 
ures" that really look like normal operation. If not excluded, such 
failures would artificially raise the missed alarm rates for all algorithms. 
Although such a circumstance would not be unfair to any particular algo- 
rithm, it would tend to interfere with the interpretation of the statistics 
gathered during the simulation process. 

The threshold value should be chosen low enough so as not to exclude 
soft failures, since they are the hardest type of failure to detect and iso- 
late, and determining the performance of the competing algorithms in 
their presence is essential. 

It can be shown that, if the standard deviation of the original gaus- 
sian distribution is <r^ and the threshold level is a, then the standard devi- 
ation of the modified distribution after the rejection process is given by 
0^2 where 

2 2 exp(-a 2 /2^) 

= exp ( -x 2 /2<t 2 )dx ' ’ 

Mode 2: Drift Rate Ramp 

The drift rate ramp is a ramp function in angular velocity. Either 
the algebraic value of the shift may be input directly or it may be Monte 
Carloed, If it is Monte Carloed, a standard deviation and threshold value 
are input. 
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Mode 3: Random Drift 


The random drift is modeled as white noise in angular velocity. 

The random drift failures are modeled the same way with a larger stan- 
dard deviation. Either the standard deviation may be input directly or it 
may be Monte Carloed. If it is Monte Carloed, a standard deviation and 
threshold value are input. 

Mode 4: Scale Factor Shift 

The scale factor shift is a step function. Its parameters are simi- 
lar to those of Modes 1 and 2. 

Mode 5: Mass Shift 

A mass shift is modeled as a change in acceleration sensitivity about 
the input and spin axes. Either the algebraic values of the two components 
of the acceleration sensitivity may be input directly or they may be Monte 
Carloed as before. If they are Monte Carloed, the same standard devia- 
tion will be used for each axis, giving circular symmetry to the error 
probability function. A threshold value is input as before, except that it 
is applied to the magnitude of the 2-vector formed by the two components 
so as to preserve the circular symmetry. 

5.1.2 Substitutional Failure Modes 

Mode 6; Zero Output 

The gyro output will be zero. 

Mode 7: Maximum Output 

Either the algebraic value of the maximum output may be input 
directly or it may be Monte Carloed. If it is Monte Carloed, the magni- 
tude of the maximum value is assigned a plus or minus sign by means 
of a uniform pseudorandom number. 

5. 2 ACCELEROMETER FAILURES 

5. 2. 1 Additive Failure Modes 

Mode 1: Bias Shift 

The accelerometer bias shift is a step function in acceleration. Its 
parameters are similar to those of gyro Mode 1. 
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Mode 2: Ramp Error 

The accelerometer ramp error is a ramp function in acceleration. 
Its parameters are simular to those of gyro Mode 2. 

Mode 3: Random Error 

The random error is modeled as white noise in acceleration. Its 
parameters are similar to those of gyro Mode 3. 

Mode 4: Scale Factor Shift 

The scale factor shift is a step function. Its parameters are simi- 
lar to those of gyro Mode 4. 

5.2.2 Substitutional Failure Modes 

Mode 6; Zero Output 

The accelerometer output will be zero. 

Mode 7: Maximum Output 

The maximum output mode is similar to gyro Mode 7. 
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6. FAILS IM 


The simulation computer program, FAILSIM, is written in 
FORTRAN IV for the GDC 6400/6500 computer. It was prepared by means 
of the TRW Timeshare System (TRW/TSS). 

6. 1 PROGRAM STRUCTURE 

Figure 6-1 shows the FAILSIM flow diagram. At the beginning of 
a run, the data are preset. Then the data for the first case are read in. 
When no data are found, the run ends. Eight separate subroutines are used 
to calculate those quantities in each FDIC algorithm that remains constant 
throughout a case. The remainder of the program is then initialized for 
the next case. 

The outer loop is on Monte Carlo trials. Nine separate subroutines 
are used to initialize the FDIC algorithms and strapdown algorithms. The 
remainder of the program is then initialized for the next trial. 

The next loop is on major cycles. It counts the number of major 
cycles in the trial. The next loop is on trajectory steps. The trajectory 
steps are the periods over which no trajectory discontinuities are per- 
missible. The innermost loop is on minor cycles. (For example, on a 
typical case the major cycle time is 1 sec, and the minor cycle time is 
1/8 sec. The trajectory was obtained from a tape having data points 
spaced 1/2 sec, so that discontinuities can occur only at the 1 /2-sec 
points. The minor cycle loop counts four minor cycles per trajectory 
step, the trajectory step loop counts two steps per major cycle, and the 
major cycle loop counts 400 major cycles per 400-sec flight.) 

The failure programmer generates the instrument errors caused by 
failures. The trajectory generator generates the nominal acceleration 
and angular velocity and combines them with the vibrational acceleration 
and angular velocity. The unfailed instrument errors are generated. 
Fourth-order Runge-Kutta-Gill integration is used to obtain the ideal 
attitude and inertial velocity and the unquantized instrument outputs. The 
sampling, quantization, and substitutional failures are then processed. 

Nine separate subroutines are used to perform the minor cycle calculations 


102 




Figure 6-1. FAILSIM Flow Diagram 
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which include the different FDIC algorithms and identical strapdown 
minor cycle algorithms. 

Nine separate subroutines are used to perform the strapdown major 
cycle algorithms. The attitude and inertial velocity errors are written 
on the plot tape. 

At the end of a trial, its statistics are stored. The trial results 
are printed out. 

At the end of a case, the statistics are processed and printed out. 
The data for the next case are then read in. 

The program has been written with a modular structure so that any 
algorithm can be readily removed for use elsewhere. Each algorithm 
comprises four subroutines: 

CONSTi: evaluate constants 

INITi: initialize variables 

MAJi: major cycle calculations 

MINi: minor cycle calculations 

(There are two exceptions: the nominal algorithm has no constants so that 
there is no CONSTO, and the Minimax Algorithm has a fifth subroutine, 
SPHERE, that performs the Bounding Sphere Algorithm when it has been 
selected. ) 

In each of the minor cycle subroutines MINi, the gyro and accel- 
erometer FDICs are coded separately to avoid the loss of time that 
would have occurred in the calling sequence if a single subroutine were 
used for both. The subroutine SPHERE is an exception to this policy. 

In addition to the algorithm subroutines described above, FAILSIM 
is divided up into subprograms which will now be described. 

FAIL . FAIL is the main program. It embodies the structure de- 
scribed above and all of the computations not relegated to the other 
subprograms. 

DERIV 1 ♦ DERIY1 is a subroutine that calculates the derivatives of 
the ideal attitude parameters (quaternions or Euler parameters) and the 
ideal inertial velocity. 
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DERIV2. DERIV2 is a subroutine that calculates the derivatives of 
the instrument outputs. The package inputs are resolved into instrument 
axes. The ideal instrument inputs are found and modified by both unfailed 
errors and additive failures to obtain the derivatives of the actual outputs. 

EVENT. EVENT is a subroutine that notes changes in the instru- 
ment states, prints out messages when changes occur, and accumulates 
statistics on the changes. 

FKB. FKB is a function that solves for (Paragraph 4. 4. 2) 

GAUSS . GAUSS is a function that converts uniform pseudorandom 
numbers into gaussian pseudorandom numbers. 

ORIENT. ORIENT is a subroutine that computes the matrices A, C, 
and F of Paragraph 3. 1. 2, the matrix R, and the matrices AG1 to AG6 
and AA1 to AA6 of Subsection 4. 1 and their products with R, and the masks 
AMASK and OMASK used in the various algorithms. 

PACKAGE. PACKAGE is a subroutine that reads in the external 
strapdown package instrument output tape data and unbuffers it. If the 
external package has three accelerometers, six pseudo accelerometer 
outputs are generated for compatibility with FAILSIM. 

PLOT. PLOT is a subroutine that finds the orientation errors in 
arc sec and the inertial velocity errors in m/sec and writes them out on 
the plot tape. 

QUAT. QUAT is a subroutine that calculates a quaternion for a 
given angle and coordinate axis and premultiplies it by an input quaternion 
to obtain an output quaternion. 

READIN . READIN is a subroutine that reads in the program inputs, 
converts units for inputs not directly involving the algorithms themselves, 
and prints out all of the inputs. (Algorithm inputs have their units con- 
verted in the CONSTi subroutines. ) 

t SYMINV . SYMINV is a subroutine that inverts a symmetric matrix 
of order 3 (only the upper triangular elements of the input matrix are 
required; all elements of the output matrix are supplied). 

TMAT. TMAT is a subroutine that calculates a rotation matrix for 
a given angle and coordinate axis and post -multiplies it by an input matrix 
to obtain an output matrix. 
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TRA J . TRAJ is a subroutine that calculates the nominal angular 
velocity and nominal acceleration for a specific trajectory. 

6. 2 INPUTS 

With the exception of taped data from an external strapdown package 
(Paragraph 6. 2. 6) and the nominal trajectory subroutine (Paragraph 6. 2. 8), 
all of the program inputs are entered by means of the NAMELIST feature. 
At the beginning of a run, all of the data is preset to 0 except IRDM, 
which is preset to 1. None of the input data is changed during a case, 
with the exception of IRDM, which always indicates the next pseudorandom 
number. 

Only those inputs which differ from 0 (or 1 in the case of IRDM) 
need be read in for the first case of a run. For subsequent cases in the 
run, only those inputs that differ from the inputs of the preceding case 
need be read in, with the exception of IRDM. If IRDM is desired to be 
the same for all cases, it must be included in the inputs to each case. 
Otherwise, each case will begin with a new pseudorandom number. In 
case it is desired to begin a run with the pseudorandom number from the 
end of a preceding run, IRDM is printed out at the end of each trial. A 
typical set of data is shown in Appendix A. 

6. 2. 1 Program Control (Subsection 6. 1) 

Variable Units Description 


DT sec 

Minor cycle period 

IRDM 

Initial random number 
(odd integer) 

LMONTE 

0 = deterministic 
failures 


1 = Monte Carlo 
failures 

LPLOT 

0 = no plot tape made 

1 = plot tape made 

N1MAX 

Number of Monte Carlo 
trials per case 

N2MAX 

Number of major cycles 
per trial 
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Variable 


Units 


Description 


N3MAX Number of trajectory- 

steps per major cycle 

N4MAX Number of minor cycles 

per step 

6.2.2 Vehicle Vibrational Motions (Paragraph 4.4. 2) 


Variable 


Units 

Description 

AVBS(ls3 ) 

rad/ 

2 

sec 

Rotational acceleration, 
vibration, body, standard 
deviation 

AVBS(4:6) 

, 2 
m/ sec 

Translational acceleration 
vibration, body, standard 
deviation 

TAUV (1:6) 

sec 


Vibration time constants 

KV(1:6) 

sec 

1 

Vibration velocity 
feedback gain 

KP(1:6) 

sec 

2 

Vibration position 
feedback gain 

RADIUS 

m 


Distance from vehicle 
c. m. to strapdown 
package 

Instrument Orientation (Subsection 4. 1) 


Variable 


Units 

Description 

IR1 



Axis of first rotation for 
R matrix 

THTR1I 

deg 


Angle of first rotation for 
R matrix 

IR2 



Axis of second rotation 
for R matrix 

THTR2I 

deg 


Angle of second rotation 
for R matrix 

IR3 



Axis of third rotation for 
R matrix 

THTR3I 

deg 


Angle of third rotation for 
R matrix 
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Variable 


Units 


THTAI(1:6) 

deg 

THTGI(ls6) 

deg 

Unfailed Errors 

(Subsection 4.2) 

Variable 

Units 

SIGUFEA(l) 

V-g 

SIGUFEA(2) 

( H-g ) ^ / (rad/sec) 

SIGUFEA(3) 

ppm 

SIGUFEA(4) 

arc sec 

SIGUFEA(5) 

arc sec 

SIGUFEG(l) 

deg/hr 

SIGUFEG(2) 

(deg/hr) 2 / 

(rad/sec) 

SIGUFEG(3) 

ppm 

SIGUFEG(4) 

arc sec 

SIGUFEG(5) 

arc sec 

SIGUFEG(6) 

deg/hr/g 

SIGUFEG(7) 

deg/hr/g 

SIGUFEG(8) 

deg/hr/g 

SIGUFEG(9) 

2 

deg/hr/g 

QACCI 

cm/ sec/pulse 

QGYROI 

arc sec/pulse 
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Output axis location angles 
for accelerometers 


Output axis location angles 
for gyros 



Misalignment 

Misalignment 

Bias 


Random 


Scale factor 
Misalignment 
Misalignment 
Acceleration sensitivity 
Acceleration sensitivity 
Acceleration sensitivity 
Anisoelastic drift 
Accelerometer quantization 
Gyro quantization 



6.2.5 Failures (Section 5) 



Up to 14 failures or healings are 

possible in a case. 

Variable 


Description 

NFAIL(1:14) 


Number of minor cycles 
before failure or healing 

NINST(1:14) 


Number of failed or 
healed instrument 

NINST = 1 j 6 failure of gyro No. NINST 

NINST = 7:12 failure of accelerometer No. NINST-6 

M0DE(1:14) 


Failure mode 

For gyros MODE = 

0 

Healing 


1 

Bias drift rate shift 


2 

Drift rate ramp 


3 

Random drift 


4 

Scale factor shift 


5 

Mass shift 


6 

Zero output 


7 

Maximum output 

For acc. MODE = 

0 

Healing 


1 

Bias shift 


2 

Ramp error 


3 

Random error 


4 

Scale factor shift 


5 

— 


6 

Zero output 


7 

Maximum output 

FAIL1I(1:14) 


First failure parameter 

FAIL2I(1:14) 
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Second failure parameter 



The units and definitions of FAIL1I and FAIL2I depend upon the 
corresponding values of NINST and MODE, and upon the value of 
LMONTE. For LMONTE = 0 (deterministic failures) and NINST = 1:6 
(gyros) 


MODE = 0, 

FAIL1I 

Not required 


FAIL2I 

Not required 

MODE = 1, 

FAIL 11 

deg /hr, value of bias 



drift rate shift 


FAIL2I 

Not required 

MODE = 2, 

FAIL1I 

2 

deg /hr , value of drift 



rate ramp 


FAIL2I 

Not required 

MODE = 3, 

FAIL1I 

2 

(deg/hr) /(rad/sec), 



random failure p. s. d. 


FAIL2I 

Not required 

MODE = 4, 

FAIL 11 

ppm, value of scale factor 



shift 


FAIL2I 

Not required 

MODE = 5, 

FAIL1I 

deg/hr/g, mass shift 



along spin axis 


FAIL2I 

deg/hr/g, mass shift 



along input axis 

MODE = 6, 

FAIL 11 

Not required 


FAIL2I 

Not required 

MODE = 7, 

FAIL1I 

deg/sec, value of 



maximum output 


FAIL2I 

Not required 

LMONTE = 0 (deterministic failures) and NINST = 7:12 (acc. ) 

MODE = 0, 

FAIL1I 

Not required 


FAIL2I 

Not required 

MODE = 1, 

FAIL1I 

pg, value of bias shift 


FAIL2I 

Not required 

MODE = 2, 

FAIL1I 

pg/hr, value of ramp 



error 


FAIL2I 

Not required 

MODE = 3, 

FAIL1I 

2 

(pg) /(rad/sec), random 



failure p. s, d. 


FAIL2I 

Not required 
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MODE = 4, 

FAIL1I 

ppm, value of scale factor 



shift 


FAIL2I 

Not required 

MODE = 6, 

FAIL1I 

Not required 


FAIL2I 

Not required 

MODE = 7, 

FAIL1I 

g, value of maximum 



output 


FAIL2I 

Not required 

LMONTE = 1 (Monte Carlo failure) and NINST = 1:6 (gyros) 

MODE = 0, 

FAIL1I 

Not required 


FAIL2I 

Not required 

MODE = 1, 

FAIL1I 

deg/hr, s.d. of bias drift 



rate shift 


FAIL2I 

deg/hr, threshold 

MODE = 2, 

FAIL1I 

2 

deg/hr , s. d. of drift 



rate ramp 


FAIL2I 

deg/hr^, threshold 

MODE = 3, 

FAIL1I 

2 

(deg/hr) /(rad/sec) s.d. 



of random failure p. s. d. 


FAIL2I 

(deg/hr)2/(rad/ sec), 



threshold 

MODE = 4, 

FAIL1I 

ppm, s.d. of scale factor 



shift 


FAIL2I 

ppm, threshold 

MODE = 5, 

FAIL 11 

deg/hr/g, s.d. of mass 



shift 


FAIL2I 

deg/hr/g, threshold 

MODE = 6, 

FAIL 11 

Not required 


FAIL2I 

Not required 

MODE = 7, 

FAIL1I 

deg/sec, magnitude of 



maximum output 


FAIL2I 

Not required 

LMONTE - 1 (Monte Carlo failures) and NINST = 7:12 (acc. ) 

MODE = 0, 

FAIL1I 

Not required 


FAIL2I 

Not required 

MODE = 1, 

FAIL1I 

pg, s.d. of bias shift 


FAIL2I 

pg, threshold 


ill 



MODE = 2, 

FAIL1I 

pg/hr, s.d. of ramp 
error 


FAIL2I 

pg/hr, threshold 

MODE = 3, 

FAIL1I 

(H L g) t */(rad/sec), s.d. of 
random failure p. s. d. 


FAIL2I 

(pg)^/(rad/sec), threshold 

MODE = 4, 

FAIL1I 

ppm, s.d. of scale factor 
shift 


FAIL2I 

ppm, threshold 

MODE = 6, 

FAIL 11 

Not required 


FAIL2I 

Not required 

MODE = 7, 

FAIL1I 

g, magnitude of maximum 
output 


FAIL2I 

Not required 

External Strapdown 

Package (Section 4. 

11 

Variable 

Units 

Description 

LPACK 


0 = Internal strapdown 


package 

1 = External strapdown 

package, 6 gyros and 
3 accelerometers 

2 = External strapdown 

package, 6 gyros and 
6 accelerometers 

(The following inputs are not required if LPACK is zero). 


IRHOl 


Axis of first rotation for 
initial quaternion 

THTRH1I 

deg 

Angle of first rotation for 
initial quaternion 

IRHQ2 


Axis of second rotation for 
initial quaternion 

THTRH2I 

deg 

Angle of second rotation 
for initial quaternion 

IRH03 


Axis of third rotation for 
initial quaternion 

THTRH3I 

deg 

Angle of third rotation for 
initial quaternion 
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Variable 

Units 

Description 

LATI 

deg 

Astronomic latitude of 
test table 

G 

/ 2 
m/sec 

Acceleration of gravity at 
test table 

OMEGA 

rad/ sec 

Earth rate 


The compensated instrument outputs in units of rad and m/sec must 
appear on a binary format magnetic tape to be read by logical unit 15. 
Each record on the tape must contain 505 or fewer words. The first 
word is an integer less than or equal to 504 and equal to the number of 
words following it in the record. If LPACK equals 1, there may be up 
to 56 nine-word data sets in the record, each data set containing the six 
gyro outputs followed by the three accelerometer outputs. If LPACK 
equals 2, there may be up to 42 twelve-word data sets in the record, each 
data set containing the six gyro outputs followed by the six accelerometer 
outputs . 

6.2.7 Algorithms (Subsection 3. 2) 

Nominal Algorithm 


Variable 

LALGO 

Adaptive 66 Algorithm 
Variable 
LALG1 


0 = 
1 = 

Units 

0 = 
I = 


Description 

do not use algorithm 
use algorithm 

Description 

do not use algorithm 
use algorithm 


113 



(The following inputs are not required if LALG1 is zero) 


EPSOA 

H-g 

Unfailed acc. error s.d. 

TT'T^O O /-« 

ii/rou 

J /I. ... 

ucg / ixr 

Unfailed gyro error s.d. 

TAUPAQ(l) 

sec 

Acc. prefilter time 
constant 

TAUPGQ(l) 

sec 

Gyro prefilter time 
constant 

CRIT1 


Loop termination 
criterion 

Fifteen Threshold Algorithm 


Variable 

Units 

Description 

LA LG 2 


0 = Do not use algorithm 

1 = Use algorithm 

following inputs 

are not required 

if LALG2 is zero) 

THAI 

Kg 

Acc. threshold 

THGI 

deg/hr 

Gyro threshold 

TAUPAQ(2) 

sec 

Acc. prefilter time 
constant 

TAUPGQ(2) 

sec 

Gyro prefilter time 
constant 

Squared Error 

Algorithm 


Variable 

Units 

Description 

LALG3 


0 = Do not use algorithm 

1 = Use algorithm 

following inputs are not required 

if LALG3 is zero) 

AK1I 

H-g 

Acc. TSE threshold, 
first failure 

AK2I 


Acc. SE/TSE threshold, 
first failure 

AK3I 

Kg 

Acc. TSE threshold, 
second failure 
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Variable 

Units 

Description 

AK4I 


Acc. SE/TSE threshold, 
second failure 

GK1I 

deg/hr 

Gyro TSE threshold, first 
failure 

GK2I 


Gyro SE/TSE threshold, 
first failure 

GK3I 

deg/hr 

Gyro TSE threshold, 
second failure 

GK4I 


Gyro SE/TSE threshold, 
second failure 

TAUPAQ(3) 

sec 

Acc. prefilter time 
constant 

TAUPGQ(3) 

sec 

Gyro prefilter time 
constant 

AKGDI 

(deg/sec) * 

Gyro dynamic error 
compensation gain 

LHEAL 


0 = no reinstatement 
permitted 



1 = reinstatement 
permitted 

Bayesian Decision Theory Algorithm 


Variable 

Units 

Description 

JLALG4 


0 = do not use algorithm 

1 = use algorithm 

following inputs 

are not required if LALG4 is zero). 

THA1I 


Acc. threshold, first 
failure 

THA2I 

f^g 

Acc. threshold, second 
failure 

THG1I 

deg/hr 

Gyro threshold, first 
failure 

THG2I 

deg /hr 

Gyro threshold, second 
failure 
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Variable 


Units 


Description 


TAUPAQ(4) 

sec 

Acc. prefilter time 
constant 

TAUPGQ(4) 

sec 

Gyro prefilter time 
constant 

Maximum Likelihood Algorithm 


Variable 

Units 

Description 

LALG5 


0 = do not algorithm 

1 = use algorithm 

SIGQAO 

Kg 

Unfailed acc. error s.d. 

SIGQGO 

deg/hr 

Unfailed gyro error s.d. 

TAUPAQ(5) 

sec 

Acc. pre filter time 
constant 

TAUPGQ( 5) 

sec 

Gyro Prefilter time 
constant 

CRIT5 

Minimax Algorithm 


Loop termination criterion 

Variable 

Units 

Description 

LALG6 


0 = do not use algorithm 

1 = use algorithm 

following inputs are not required if LALG6 is zero). 

LSQA 


0 = Acc. use Bounding 

Sphere Algorithm 

1 = Acc. use Least- 

Squares Algoritnm 

LSQG 


0 = Gyros use Bounding 

Sphere Algorithm 

1 = Gyros use Least- 

Squares Algorithm 

THAMI 

Kg 

Max. filtered error of 
unfailed acc. 

THASI 

Kg 

Max. unfiltered error of 
unfailed acc. 
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Variable 

Units 

Description 

THGMI 

deg/hr 

Max. filtered error of 
unfailed gyro 

THGSI . 

deg/hr 

Max. unfiltered error of 
unfailed gyro 

TAUPAQ(6) 

sec 

Acc. prefilter time 
constant 

TAUPGQ(6) 

sec 

Gyro prefilter time 
constant 

(THASI is not required if LSQA equals 1 
equals 1). 

; THGSI is not required if LSQG 

Sequential Algorithm 


Variable 

Units 

Description 

LALG7 


0 = do not use algorithm 

1 = use algorithm 

(The following inputs 

are not required if LALG7 is zero). 

SIGPA(l) 

^g 

Acc. input state initial 
s. d. 

SIGFA(2) 

^g 

Acc. error state initial 
s. d. 

SIGQA(l) 

^g 

Acc. input state noise 
s. d. 

SIGQA(2) 

Hg 

Acc. error state noise 
s. d. 

SIGRA 


Acc. measurement error 
s. d. 

TAUSA(l) 

sec 

Acc. input s.t.m. time 
constant 

TAUSA(2) 

sec 

Acc. error s.t.m. time 
constant 

THRAI(l) 

l^g 

Acc. high level threshold 

THRAI(2) 


Acc. low level threshold 
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Variable 

Units 

Description 

SIGFG(l) 

deg/hr 

Gyro input state initial 
s. d. 

SIGPG(2) 

deg/hr 

Gyro error state initial 
s. d. 

SIGQG(l) 

• deg /hr 

Gyro input state noise 
s. d. 

SIGQG(2) 

deg/hr 

Gyro error state noise 
s. d. 

SIGRG 

deg /hr 

Gyro measurement error 
s. d. 

TAUSG(l) 

sec 

Gyro input s. t. m. time 
constant 

TAUSG(2) 

sec 

Gyro error s.t. m. time 
constant 

THRGI(l) 

deg/hr 

Gyro high level threshold 

THRGI(2) 


Gyro low level threshold 


6. 2. 8 Nominal Trajectory (Paragraph 4. 4. 1 and Appendix A) 

The nominal trajectory is specified by a subroutine named TRAJ 
written in FORTRAN IV. A shuttle boost trajectory is included in the 
program at present. The trajectory is broken up into piecewise segments 
of analytical functions of time. The logic to choose the correct segment 
must be supplied. It may depend upon the time T, provided there is no 
discontinuity in the functions at the switching time. Any values of T are 
permissible for use in the segment selection logic. If there are discon- 
tinuities present, they may only occur at the times 

tj. = N • N4MAX • DT N = 0,1,2,... (6-1) 

The number of the discontinuity opportunity is given by N 

N = N3MAX '• (N2 - 1) + N3 (6-2) 

where N2 is the major cycle counter and N3 is the trajectory step counter. 
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Eq. 6-2 should be the first executable statement of the subroutine. The 
logic to choose the correct segment must depend upon N if there is a 
discontinuity. The value of T corresponding to N is given by Eq. 6-1. 
Logic based on T may be mixed with logic based on N. The body angular 
velocity is specified by VPAB(1:3), and the body acceleration is specified 
by APLB(1:3). It is convenient to take the body axes as: 

(1) Roll (forward) 

(2) Pitch (right) 

(3) Yaw (down). 

6.3 OUTPUTS 

A typical run is shown in Appendix A. The outputs will be de- 
scribed here. 

6, 3. 1 Inputs 

At the beginning of the output we find the program name, the date, 
and the time. Then all of the NAMELIST inputs are printed out in groups. 
If LALGi is 0, the remaining inputs for algorithm i are not printed, even 
if they were present in the input data. 

6. 3. 2 Events 

Following the inputs, we find the output data for each trial. First 
we see a listing of all of the events taking place during the trial. When- 
ever an instrument is failed, the time, instrument type, instrument num- 
ber, and mode of the failure are given. Whenever an instrument is 
healed, the time, instrument type, and instrument number are given. 
Whenever an algorithm changes state, the time, algorithm number, 
instrument type, instrument number, and type of state change are given. 
There are nine possible types of state change of interest: 

ALL CLEAR The algorithm changes the state of 

a healed instrument from failed to 
unfailed 

CORRECTED FALSE The algorithm changes the state of 

ALARM an unfailed instrument from failed 

to unfailed 

FALSE ALARM The algorithm changes the state of 

an unfailed instrument from unfailed 
to failed 
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LOST ALARM 

The algorithm changes the state of 
a failed instrument from failed to 
unfailed 

TRUE ALARM 

The algorithm chanuea the state of 

o o 

a failed instrument from unfailed 
to failed 

FORTUITOUS TRUE 
ALARM 

An instrument fails when the 
algorithm has already determined 
its state to be failed 

FORTUITOUS ALL 
CLEAR 

An instrument is healed when the 
algorithmhas already determined 
its state to be unfailed 

MISSED ALARM 

The trial terminates and the 
algorithm has determined the state 
of a failed instrument to be unfailed 

MISSED ALL CLEAR 

The trial terminates and the 
algorithmhas determined the state 
of a healed instrument to be failed. 


Figure 6-2 shows the state transitions and the corresponding messages. 
The left-hand two states (circles) correspond to an unfailed instrument, 
the middle two to a failed instrument, and the right-hand three to a healed 
instrument. The upper four states correspond to an instrument state in 
the algorithm of failed (alarm), and the lower three states correspond to 
a state of unfailed (no alarm). Horizontal transitions are caused by fail- 
ure or healing, curved line transitions by changes of algorithm instru- 
ment state, and vertical straight line transitions by the end of the trial. 
The upper right-hand state is split so as to distinguish between all clears 
and corrected false alarms. 

Unlike the other algorithms, the Adaptive 66 and Maximum Likeli- 
hood Algorithms do not have internal logical states that give obvious 
meanings to the concept of alarm or no alarm. Rather than obtain no 
such information for these algorithms, it was decided to invent an arbit- 
rary criterion for the presence or absence of an alarm. This criterion, 
for both algorithms, is that an alarm occurs when the variance of an 
instrument is found to be greater than three times the variance of unfailed 
instruments. 
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END OF TRIAL 



Figure 6-2. Event Logic 

6.3.3 Figures of Merit 

The figures of merit are printed out at the end of the trial. They 
appear in nine columns for algorithms 0 to 8, and four rows. The first 
row gives the magnitudes of the small angle error vectors between the 
attitudes of algorithms 0 to 8, and the ideal attitude. The second row 
gives the magnitudes of the small angle error vectors between the atti- 
tudes of algorithms 1 to 8, and the attitude of algorithm 0, the nominal 
algorithm. The third row gives the magnitudes of the error vectors 
between the inertial velocities of algorithms 0 to 8, and the ideal inertial 
velocity. The fourth row gives the magnitudes of the error vectors 
between the inertial velocities of algorithms 1 to 8, and the inertial 
velocity of algorithm 0. The first two rows are in arc sec and the last 
two in m/ sec. 

The second and fourth rows are useful in discerning differences 
between the algorithms when their answers are close to the nominal 
algorithm, that is when the errors caused by failures are being obscured 
by the effects of strapdown algorithm truncation error and unfailed errors. 
The absolute significance of these figures of merit is less clear. 


121 



The algorithm numbering scheme is: 


0 Nominal 

1 Adaptive 66 

2 Fifteen Threshold 

3 Squared Error 

4 Bayesian Decision Theory 

5 Maximum Likelihood 

6 Minimax 

7 Sequential 

8 Not yet implemented 
(intended for Adaptive 72). 

6.3.4 Algorithm Execution Times 

Following the figure of merit, we find the total times (in sec) spent 
by the central processor in the FDIC portion of each algorithm during the 
trial, and the peak times required by each algorithm in any single minor 
cycle. Unfortunately, the latter are quantized too coarsely (msec) to be 
of much value. 

The final random number and the total central processor time used 
so far in the trial are then printed out. 

6.3.5 Statistics 

At the end of the case, the statistics of the trials are printed out. 
The mean magnitudes of the figures of merit are printed out in the same 
format as described above. The total FDIC time for all trials and the 
peak FDIC times for any trial are also printed out. A summary of all 
the events occurring during the case is also presented. 

6. 3. 6 Plot Tape 

The plot tape is written in binary format on logical unit 13 using 

binary blocking. If binary blocking is not desired, or if the program is 

to be used with other than CDC 6X00 computers, a DATA statement and 

the third executable statement in the main program, FAIL, should be 

removed. These statements are: 

DATA IRAY/13/ 

CALL FTNBIN(1, 1, IRAY) 
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Each record on the plot tape contains 19 words. The first word is the 

til 

time. The (2i + 2) word is the magnitude of the small angle attitude 

til 

error vector in arc sec for algorithm i; the (2i + 3) 1 word is the mag- 
nitude of the inertial velocity error vector in m/sec for algorithm i; 
i = 0, 8. 
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7. SETTING THE PARAMETERS 


The seven algorithms actually implemented have among them 
sons 26 adjustable parameters for the gyrus etiiu. ion cne acceiero— 
meters; the Sequential Algorithm accounts for a third of these. Setting 
these parameters correctly presents quite a problem. Ideally, the 
parameters should be set optimally for each algorithm. However, even 
reaching an acceptable definition of optimality is a difficult task. 

The analytical setting of the parameters presents formidable dif- 
ficulties, since all of the algorithms are nonlinear, and a different analy- 
sis would have to be performed for each. In fact, if one were capable of 
setting the parameters analytically, this would imply the ability to deter- 
mine performance analytically, and FAILS IM would be unnecessary. 

One example of the type of difficulty encountered in trying to predict 
behavior analytically is the "failure induced false alarm. " A configuration 
of unfailed errors that could not produce a false alarm by themselves 
may, in the presence of a failure of one instrument, cause a false alarm 
for a different instrument. If an instrument is already failed and isolated, 
the false alarm precludes isolation of the new failure. 

Table 7-1 shows the number of parameters in each algorithm. It 
also shows whether or not independent threshold settings are available for 
the second failure. 


Table 7-1. Adjustable Parameters in the Several Algorithms 


Algorithms 

Parameters 

Threshold for 
Second Failure 

Gyros 

Accs . 

Adaptive 66 

2 

2 

Same as first 

Fifteen Threshold 

2 

2 

Same as first 

Squared Error 

6 

5 

Independent 

Bayesian Decision Theory 

3 

3 

Independent 

Maximum Likelihood 

2 

2 

Same as first 

Minimax/ Least Squares 

2 

2 

4. 2 dB above first 

Adaptive 72 

? 

? 

Independent 

Sequential 

9 

9 

Same as first 
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Happily, the authors of the Squared Error Algorithm have 
supplied four of the thresholds: 

AK2I = GK2I = 4/9 (7-1) 

AK4I = GK4I = 8/21 (7-2) 

If we assume that setting AKGDI to zero has no major effect on perfor- 
mance, then we need determine only three parameters for the gyros and 
three for the accelerometers. Then only the Sequential Algorithm has 
more than three parameters to determine. The difficulty in selecting 
parameters for this algorithm has been alleviated somewhat by the analy- 
sis of Appendix B. By its use, one can assign the parameters of the 
Sequential Algorithm to give an equivalent steady-state time constant. 

To provide a reasonably fair way of setting the parameters, it was 
decided to allot equal amounts of effort to each algorithm. After all, if 
an algorithm requires less effort than another in parameter setting, it 
has an advantage that ought to show up in the final results . 

After consideration of various possibilities, it was decided that a 
reasonable approach, capable of being achieved in a limited number of 
computer runs, would be to set all of the thresholds for the same false- 
alarm rate. Thus all algorithms would give about the same performance 
in the absence of failures (the most probable situation for any particular 
mission), and their relative performance in the presence of failures 
would serve to establish their relative merits. 

The procedure followed was to choose two values for the time con- 
stants: 100 and 10 sec. The 100-sec value was used first. A series of 
runs was made with no failures. The thresholds were adjusted in steps 
of 2 dB until a value was found for each algorithm for which no false alarm 
alarms occurred while a 2-dB lower threshold produced at least one 
false alarm. (The accelerometer threshold-setting process was not com- 
pleted for lack of time, so the following comments apply to the gyros 
only. ) 
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Runs were then made with small MODE = 1 failures at 200 sec 
(halfway through boost). False alarms occurred in the Fifteen Threshold, 
Minimax, and Sequential Algorithms. Raising the thresholds 2 dB elimi- 
nated the false alarms. 

More runs were made with a large MODE = 1 failure at 190 sec and 
a small MODE = 1 failure at 200 sec. The purpose of these runs was to 
set the thresholds for second failures in those algorithms that had them. 
The large failure was chosen large enough to cause all of the algorithms 
to isolate the failure within one minor cycle, so that no data from the failed 
instrument were used. It was found necessary to set the second -failure 
thresholds higher than the first- failure thresholds in the Squared Error 
and Bayesian Decision Theory Algorithms in order to eliminate false 
alarms. For the same reason, it was necessary to raise the thresholds 
of the Adaptive 66, Fifteen Threshold, Maximum Likelihood, and Sequen- 
tial Algorithms, thus degrading their performance on first failures. Thus 
the desirability of independent first- and second-failure thresholds seems 
obvious. However, it was not found necessary to raise the single thresh- 
old of the Minimax Algorithm, apparently because the built-in 4. 2 dB 
ratio between effective first- and second-failure threshold levels is close 
to optimal. For the Maximum Likelihood and Sequential Algorithms, the 
thresholds had to be raised so high in order to eliminate false alarms 
that a true alarm was not obtained on the second failure. 

This procedure was repeated with a 10-sec time constant. With 
this value of time constant, the Maximum Likelihood and Sequential 
Algorithms succeeded in obtaining true alarms without false alarms. The 
10-sec time constant was found to be somewhat better for all algorithms, 
so it was chosen for the remaining runs. The thresholds selected to go 
with the 10-sec time constant were considered the final set for use in the 
competition. However, the first double -failure test run produced false 
alarms (before any failures) for the Squared Error, Bayesian Decision 
Theory, and Minimax Algorithms, so their thresholds were all raised by 
2 dB, eliminating the false alarms. 
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8. SIMULATION RESULTS 


The symbols used to plot the errors of the various algorithms are: 

O Nominal Algorithm or multiple point 
A Adaptive 66 Algorithm 
+ Fifteen Threshold Algorithm 
X Squared Error Algorithm 
<•> Bayesian Decision Theory Algorithm 
Q Maximum Likelihood Algorithm 
□ Minimax Algorithm 
Z Sequential Algorithm 

8. 1 100-SEC TIME CONSTANT RUNS 

The gyro FDIC algorithm parameters for the 100- sec time constant 


runs are given in Table 8-1. (The threshold setting process was not suc- 
cessful for the Maximum Likelihood and Sequential Algorithms so they do 
not appear in Table 8-1.) 

Table 8-1. FDIC Parameters 

Adaptive 66 


EPSOG = 1. 585 

TAUPGQ(l) = 100. 

CRIT1 = 0. 1 


Fifteen Threshold 


THG1 = 1. 585 

TAUPGQ(2) = 100. 

Squared Error 


GK1I = 1. 41 

GK2I =4/9 

GK3I = 2. 51 

GK4I = 8/21 

TAUPGQ(3) = 100. 

AKGDI = 0. 

LHEAL = 0 


Bayesian Decision Theory 


THG1I = 1. 995 

THG2I = 3. 98 

TAUPGQ(4) = 100. 


Minimax 


THGMI = 0. 617 

TAUPGQ(6) = 100. 

LSQG = 1 
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8. 1. 1 First Failures 

A series of runs was made with a MODE=l gyro failure of varying 
magnitude at 200 sec using the parameters shown above. Figure 8-1 
shows the attitude error at the end of the single trial in arc sec versus the 
magnitude of the failure in deg/hr. The same initial random number 
was used for each trial. 

8. 1. 2 Second Failures 

A series of runs was made with a large MODE=l gyro failure at 
190 sec. All of the algorithms switched out the failed gyro within one 
minor cycle, so that no erroneous data were incorporated in the strapdown 
solution. A second MODE = l gyro failure of varying magnitude occurred 
at 200 sec, using the parameters shown above. Figure 8-2 shows the 
results. 

8.2 10 -SEC TIME CONSTANT RUNS 

The gyro FDIC algorithm parameters for the 10- sec time constant 
runs are given in Table 8 -II. 

8. 2. 1 First Failures 

A series of runs was made with a MODE=l gyro failure of varying 
magnitude at 200 sec using the parameters shown above. Figure 8-3 
shows the results. 

8. 2. 2 Second Failures 

A series of runs was made with a large MODE=l gyro failure at 
190 sec and a smaller MODE = l gyro failure of varying magnitude at 
200 sec, using the parameters shown above. Figure 8-4 shows the 
results. 

The 10-sec time constant results look better than the 100-sec time 
constant results for the larger failures and somewhat worse for the 
smaller failures. The choice of 10 sec might minimize the average error, 
while the choice of 100 sec might minimize the worst-case error. An 
intermediate value of time constant might be a good compromise. In any 
event, 10 sec was chosen as the time constant for the subsequent runs. 
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Figure 8-1. Orientation Error Magnitude Versus Failure Magnitude 
First Gyro Failure, 100-Sec Time Constant 



Figure 8-2. Orientation Error Magnitude Versus Failure Magnitude 
Second Gyro Failure, 100- Sec Time Constant 
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Figure 8-3. Orientation Error Magnitude Versus Failure Magnitude 
First Gyro Failure, 10-Sec Time Constant 



Figure 8-4. Orientation Error Magnitude Versus Failure Magnitude 
Second Gyro Failure, 10-Sec Time Constant 


130 


Table 8 -II. FDIC Parameters 


Adaptive 66 


EPS0G = 3. 16 
CRIT1 = 0. 1 

TAUPGQ(l) = 10. 

Fifteen Threshold 


THG1 = 5. 01 

TAUPGQ(2) = 10. 

Squared Error 


GK1I = 3. 54 

GK2I =4/9 

GK3I = 6. 31 

GK4I = 8/21 

TAUPGQ(3) = 10. 
LHEAL = 0 

AKGDI = 0. 

Bayesian Decision Theory 


THG1I = 5. 01 
TAUPGQ(4) = 10. 

THG2I = 10. 

Maximum Likelihood 


SIGQG0 = 8. 78 
CRIT5 = 0. 1 

TAUPGQ(5) = 10. 

Minimax 


THGMI = 1. 55 
LSQG = 1 

TAUPGQ(6) = 10. 

Sequential 


SIGPG(l) = 206000 

SIGPG(2 ) = 4. 03 

SIGQG(l) = 206000 
SIGRG = 5 1 

SIGQG(2 ) = 0. 55 

TAUSG(l) = 0. 0128 

TAUSG(2) = 20. 

THRGI(l) = 618000 

THRGI(2) = 0. 196 


This behavior of the errors in response to changes in the time con- 
stant of the prefilter (with concomittant threshold changes to keep the 
false alarm rate constant) can be verified from the analysis of Subsec- 
tion 3. i qualitatively, if not quantitatively. Figure 8-5 shows Eq. 3-110 
plotted for two values of b and Tf. The values of Tf are the two values 
used in the simulation. When the time constants were decreased from 
100 to 10 sec, the thresholds were increased by 8 to 10 dB. Thus the 
ratio between the two values of b used in Figure 8-5 was taken as 9 dB. 

The left-hand portion of both curves is simply the failure magnitude multi- 
plied by 200 sec (the time from the failure to the end of boost). 
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Figure 8-5. Theoretical Error Magnitude Versus 
Failure Magnitude 


8. 3 DOUBLE SIMULTANEOUS FAILURE RUNS 

A series of runs with double simultaneous gyro failures of different 
MODE values was made as a severe test of the capabilities of the different 
algorithms. The three algorithms which gave the poorest results in the 
previous tests, the Adaptive 66, Maximum Likelihood, and Sequential Algo- 
rithms, also required the most central processor time. It was therefore 
decided to eliminate them from contention at this point, and devote the 
remaining resources to the evaluation of the more successful algorithms. 

On the first trial, the Squared Error, Bayesian Decision Theory, 
and Minimax Algorithms all experienced false alarms before the failures 
occurred. Their thresholds were all raised 2 dB, eliminating the false 
alarms. The gyro FDIC algorithms parameters for these runs are given 
in Table 8 -III. 
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Table 8-IIL FDIC Parameters 


Fifteen Threshold 


THG1 = 5. 01 

TAUPGQ(2) = 10. 

Squared Error 


GK1I = 4. 46 

GK2I = 4/9 

GK3I = 7. 94 

GK4I = 8/21 

TAUPGQ(3) = 10. 
LHEAL = 0 

AKGDI = 0. 

Bayesian Decision Theory 

• 

THG1I = 6. 31 
TAUPGQ(4) = 10. 

THG2I = 12.5 9 

Minimax 


THGMI = 1. 953 
L.SQG = 1 

TAUPGQ(6) = 10. 


Every combination of two failures of MODE 1 through MODE 6 was 
tried for a total of 21 cases. MODE 7 was omitted because it was felt that 
no algorithm would have much difficulty with so hard a failure. Different 
initial random numbers were used for each case. Only one trial was made 
per case. The magnitudes of the attitude error vector with respect to the 
ideal attitude and with respect to the nominal algorithm attitude were 
averaged over the 21 cases, giving the results in Table 8-IY. 


Table 8 -IV. Results of Double Simultaneous Failures (arc sec) 


Algorithm 

Ideal 

Nominal Algorithm 

Nominal 

434 

0 

Fifteen Threshold 

688 

396 

Squared Error 

761 

465 

Bayesian Decision Theory 

841 

563 

Minimax/Least Squares 

672 

366 
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The events occurring during the 21 trials are summarized in 
Table 8-V. The prevalence of false alarms indicates that the thresholds 
may still be set a bit too low. The prevalence of missed alarms is partly 
due to the choice of too low-level failures in some of the failure modes. 


Table 8-V. Summary of Events 


Events 

Algorithms 

Fifteen 

Threshold 

Squared 

Error 

Bayesian 

Decision 

Theory 

Minimax 

Least- 

Squares 

False alarms 

4 

5 

3 

4 

True alarms 

23 

25 

21 

25 

Missed alarms 

19 

17 

21 

17 


We see that the two algorithms using the indirect test signals 
(parity signals) do somewhat better than the two using the direct test 
signals (residuals). We also note that the test signals of the three algo- 
rithms already eliminated from contention are more closely related to the 
direct than to the indirect test signals. Thus we find a reasonably strong 
presumption in favor of the indirect test signals being the better choice. 

With the benefit of hindsight, we may find plausible reasons why 
this result is obtained. One reason is that, with 15 signals as compared 
to six, one can discriminate a greater number of different states in v- space. 
Another reason is that there is only a 7 dB ratio between the contribution 
of the i** 1 instrument error to the residual and the contribution of the 

*■» J."L. 

j instrument error to the i residual, whereas there is an infinite ratio 
between the contributions to a parity signal of the error of an instrument 
included in it as compared to an instrument excluded from it. 

8. 4 MINOR CYCLE TIME 

All of the preceding runs were made with a minor cycle time of 
12 5 msec. Originally it was intended to use 10 msec, which is a value 
more typical of practical strapdown packages. However, the running 
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time on the CDC 6400/6500 proved to be excessive. The 125-msec period 
was selected as a good compromise between running time and loss of accu- 
racy. In order to estimate the loss of accuracy, a single trial was run 
with the nominal algorithm only using the 10-msec period. The resulting 
attitude error magnitude of 113 sec may be compared with the value of 
434 sec given in Table 8-IV. 

If we RSS the 113 sec nominal algorithm error with the error with 
respect to the nominal algorithm achieved in the results above, we get a 
crude indication of the performance achievable in the presence of failures 
with the Minimax Algorithm with the 10 -msec period. 

113 2 + 366 2 = 384 2 (8-1) 

This estimate shows some 12 dB of performance degradation in the 
presence of double failures as compared to the nominal algorithm. In the 
worst case of the 21 cases, the degradation was 27 dB. Of course, such 
occurrences will probably be rare, but it is interesting to note that even 
the best algorithm can occasionally do rather poorly, (However, a 
slightly higher threshold setting might have eliminated the worst case, 
which was caused by false alarms. ) 
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9. ALGORITHM SIZING AND TIMING 
9. 1 ALGORITHM SIZING 

Algorithm sizing comparisions are made by finding the storage 
requirements for the FDIC algorithms on the CDC 6400/6500 computer. 
Only the MINi subroutines are included. The number of words used for 
instructions, the number of words used for data, and the total of the two 
are presented in Table 9-1. 


Table 9-1. Algorithm Sizing 


Algorithms 

Instruction Words 

Data Words 

Total 

Adaptive 66 

460 

156 

6 1 6 

Fifteen Threshold 

358 

120 

478 

Squared Error 

608 

153 

761 

Bayesian Decision Theory 

526 

137 

663 

Maximum Likelihood 

640 

202 

842 

Minimax 

700 

167 

867 

Sequential 

447 

360 

807 


It should be noted that CDC 6400/6500 words contain several instruc- 
tions each, so that an accurate instruction count has not been obtained. 

(The Bounding Sphere Algorithm is excluded from the Minimax Algorithm 
figures since it was not used; however, the logic deciding whether to use 
least-squares or bounding spheres and the calling sequence to the sub- 
routine SPHERE is still included. Also, the Minimax Algorithm can be 
simplified somewhat if it is used only with least-squares. Therefore, the 
figures for the Minimax Algorithm must be considered to be a bit high. ) 

9. 2 ALGORITHM TIMING 

Algorithm timing comparisons are made by keeping track of the 
central processor time spent on each algorithm. As mentioned above, 
the peak timing results are too coarsely quantized to be very useful. The 
total time used by the different algorithms in a typical run is given in 
Table 9-II. 
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Table 9-II. Total Processor Time (sec) 


Algorithms 

Time 

Adaptive 66 

10. 6 

Fifteen Threshold 

8.5 

Squared Error 

6. 3 

Bayesian Decision Theory 

7. 0 

Maximum Likelihood Theory 

15. 6 

Minimax 

14. 6 

Sequential 

188. 0 


However, the Adaptive 66 and Maximum Likelihood Algorithms, 
unlike any of the others, have iterations or loops of indefinite length in 
their structure which, on occasion, can consume much more time than 
even the Sequential Algorithm. For example, the Adaptive 66 Algorithm 
once used up 479. 3 sec of central processor time in a single trial. (The 
Maximum Likelihood Algorithm was not turned on, fortunately, or it 
might have done even worse. ) The timing on the Minimax Algorithm may 
be somewhat high, for the same reasons as given in Section 9. 1. The 
average (but not the peak) time could be decreased somewhat at the cost 
of increased size by reverting to the original scheme, which has a separate 
failure detection scheme with failure isolation not occurring unless a 
failure is detected. 
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10. AREAS FOR FUTURE STUDY 


Although the original tasks required in the study have all been 

i_j a 1 M wrl-tA-MA *s^/Kf4Ar»a1 T*rnt»lr /vnrrV»4- fn Ko rlrm/a 

COinpie l»cu, Ct lUUllUCl Ui O.J. ^«.o w x.l<s J. v. »rw**v vw 

have become apparent during the course of the study. They will now be 

discussed. 


Up to the time that the Adaptive 72 Algorithm appeared, the question 
of filtering seemed to be merely a choice between two linear low-pass fil- 
ters with slightly different impulse responses. The Adaptive 72 Algo- 
rithm uses a nonlinear filter described as "a suboptimal detection system 
based on Wald's sequential analysis using the concept of information 
value and information feedback. " This filter could be used by the Fifteen 
Threshold, Bayesian Decision Theory, and Minimax Algorithms as well 
as by the Adaptive 72 Algorithm. It would be worthwhile to compare the 
three kinds of filters in the same algorithm. 


The Adaptive 72 Algorithm failure detection, isolation, and correc- 
tion schemes should be coded for evaluation in FAILSIM when the docu- 
mentation of the SIRU application becomes available. The identification 
and recertification schemes need not be studied at this time as there is 
no competitor on the scene — merely confirming the originator's results 
would not be worthwhile. 


The original program plan did not include the evaluation of the 
effects of "glitches" (temporary failures) on algorithm performance. How- 
ever, the capability of healing failures was included in the failure pro- 
grammer of FAILSIM. A study of the effects of glitches using FAILSIM 
would be desirable. It will probably lead to the modification of some of 
the algorithms, however, because none of them appear to have been 
designed with glitches specifically in mind. The hypothesis that an instru- 
ment has healed itself is not the complement of the hypothesis that an 
instrument has failed, because, in the first case, we know which instru- 
ment we are talking about before we make the decision, while in the second 
case we do not and must, in effect, make six decisions. Thus the decision 
that an instrument should be reinstated should not necessarily be based 
on the disappearance of the conditions that caused it to be classified as 
failed in the first place. 
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A goal (rather than a requirement) of the original program plan 
was to combine the best features of the algorithms into a super -algorithm. 
In addition to the obvious combination of the best filter or prefilter with 
the best isolation technique, other possibilities exist. An example is 
that of applying the signal locking technique used with the Fifteen Thresh- 
old Algorithm to the Minimax Algorithm to see whether improved 
isolation results. Thus this goal is still a desirable one. 
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APPENDIX A 
TYPICAL RUN 

The inputs and outputs of a typical run are presented in this appendix. 
A. 1 NAMELIST INPUTS 

The NAMELIST inputs are presented in this section. They are typi- 
cal of the inputs used during the final series of runs, the double simul- 
taneous failure runs. Although algorithms 1, 5, and 7 are turned off, 
their parameters are left in the input file for convenience. The values of 
these parameters are the final set used in the 10-sec time constant runs. 

This run has only one case. To add more cases, one simply appends 
their data decks. The data for each case must be preceded by a $INPUT 
card and terminated by a $. All cards must start in column 2 , except the 
comments which have a C in column 1. 

The case has only one trial. To add more trials, set NIMAXto the 
number of trials. 
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COATA FCR FAIL IN4, CASE 1 
$ INPUT 
OT = .125 
I RDM * 251 7 56 6905 
LMCNTE = 0 
LPLQT =0 
N1MAX = 1 
N2MAX = 400 
N3MAX = 2 
N4MAX = 4 

A VBS = 1., l.t l.t .981, .981, .981 
T AUV = .04, .04, .04, .04, .04, .04 
KV = 7., 7., 7., .7, .7, .7 
KP = 25., 25., 25., .25, .25, .25 
RADIUS = 6. 

I R 1 = i, IP2 = 2, IR3 = 3 

THTR1I = 0., THTR2I = 0., THTR3I = 0. 

TFTAI = 0 . , 0., 0., O.t 0., 0. 

THTGI — 0., 0., 0., 0 » , 0., 0. 

SIGUFEA = 66., 1., 70., 50., 50. 

SIGUFEG = .063, .01, 100., 10., 10., .2, .02, .2, .04 
QACCI = 2. 

QGYROI = 5. 

NFAIL(l) = 1600, 1600 
NINST(l) = 1, 6 
MODE 111 = 1, 1 
FAILlI(l) = 10., 10. 

LALGO = 1 
LALG1 = 0 
EPSOA = 1585. 

EPSOG = 3.16 
TAUPAUU) = 10. 

TAUPGQ(l) = 10. 

CRIT1 = .1 
LALG2 = 1 
THAI = 1995. 

THGI = 5.01 
TAUPAQI2) = 1C. 

TAUPGQI2) = 1C. 

LALG3 = 1 
LHEAL = 0 


144 


AK1I = 2240. 

AK2I = • 444444444444444 
AK3 I = 1995. 

AK4I = .38 09 5238 09 52381 
GK1I = 4.46 

GK2I = .444444444444444 
GK3 I = 7.94 

GK4I = .380952380952381 

T AUPAQ( 3) * 1C. 

T AUPGQ ( 3 ) = 10. 

LALG4 = 1 
THAU = 3160. 

THA2I = 3160. 

THG1I = 6.31 
THG2I = 12.59 
T AUPAQ (4) = 1C. 
TAUPGQC4) = 1C. 

LALG5 = 0 
SIGQAO = 5540. 

SIGOGO a 8.78 
TAUPAQJ5) = 10. 

T AUPGQ ( 5 ) = 1C. 

CRIT5 = .1 
LALG6 = 1 
LSQA = 1 
LSQG = 1 
THAMI = 978. 

THASI » 100. 

THGMI = 1.953 
THGSI = .1 
T AUPAQJ6) = 10. 

T AUPGQ ( 6 ) = 1C. 

LALG7 x o 

SIGPA = 3.E6, 65.8 
S1GPG x 206000., 4.03 
SIGQA = 3.E6, 8.99 
SIGQG = 206000. t -550 
SIGRA = 833. 

SIGRG = 51. 

TAUSA = .0128, 20. 

TAUSG = .0128, 20. 

THRAI = 9.E6, 7.58 
THRGI = 618000., .196 
SEND 
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A. 2 TRAJECTORY SUBROUTINE 


The subroutine, TRAJ, that defines the nominal trajectory is pre- 
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graph 6. 2. 8 for a detailed description of the method. Note that VPAB and 
APLB are zeroed out in the main program, FAIL, before each trial, so 
that statement 80, for example, computes only the nonzero component of 
angular velocity. 
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♦ TRAJ 


SUBROUTINE TRAJ 

COMMON ROM, A(6,3), AMASK 115), BAC<3,6,5), BGQ<3,6,9), C<3,6), 

1 CF 1 15,3 ) , MPSSMCG , OMASK (6 ) , OUTPI12), PHI STRO 13,9) , RAOSEC, 

2 RH0C(4,91 .“ftHOO m r S'ACm , SGOm ,“7TTTKFre,2,2» i VBCT 3', 5) , 

3 VIQ(3,9> 

COMMON ACC (3 ) , AKGO I, AK1I, AK2I, AK3I, AKA 1 , AR{12,3»3), OT, 

T^cTMTQTTrrzrTpyoTrv t nnrrrTn ? it ta il tt tdt ~ cki t , gk7t — G or 7 

2 GKAI, I, IDATA, INST C 12 ) , IROM, IR1, IR2, IR3, K0UNTA(8,9), 

3 KCUNTG (8,9 ) , LALG(S), LHEAL, L MONTE, LPLCT , ISGA, LSGG, MGDEI1A), 

| 4 NFXnTRT* NTN5TTTTI , NIMaX, N2, N2WAX, N3, N3HAX, K^HSTT UHGT3T 

5 COT, CACC , QGYRO, CRITl, CRIT5, SAQG (9) , SAQ1 (9) , SGC0I9), 

I 6 SGQ1(9), S1GPA (2) , SIGPG(2), SIGQAI2), SIGCAO, SIGQGI2), SIGQGC, 
TsrGRA.T ^ ), “STGliFEG ( 91 , TAUPAQ (6 IT TAWGG'I 61 , “ 

8 T AUSA (2 ) , TAUSG ( 2 ) , TFAI L( 12 ) , THAI, THAMI, THASI, THAU, THA2I, 

9 THGI , THGMI, THGSI , THG1I, THG2I, THRA 1(2), THRGK2), THTA (6) , 

I A~TKT0T5TT~THTKT7 TTi T R ' 2 , ~T H T R 3 UFE 1 1 2 ,9)T Tt T ) , ' XFA1LH 14) , 

! 8 XFAIL2IU) 

COMMON APLB (3 ) , AVB (6 ) , A VB A 1 6 ) , AVBBI6), AVBS(fc), AVB016), 
I~AV0TT8TT AVB 2 ( OYAVBTT6 T,~ATr6TT~AZ 06T , “BCT 63 , TIT6T, ET2T6T, 

2 00(6), C2 ( 6) , OX (7 ) , KP(6), KV(6), 0(7), RAC10S, TAUV16), 

3 TMIA.A), VP AB ( 3 ) , VVA8(3) 

[CHPITO'CTTRMFC TRAJECTORY' FT FOR F ACC ITTnOF' RANDOM VIBRATIONS 

I N = N3MAX* (N2 - 1) ♦ N3 

C BODY ANGULAR VELOCITY 



IF 

Tn 

• GT. 3*131 

GCrTCTAC 


IF 

IN 

.GT. 22 A ) 

GC TO 20 


IF 

(N 

.GT. AO) 

GO TO 10 


IF 

TAT 

• G T • 16 1 

so, ec ______ 

Il 0 

IF 

(N 

-GT. 1A0) 

lie, ICC 

>20 

IF 

IN 

-GT. 232) 

GO TG 30 


IF 

IN 

.GT . 2251 

13<r,I2C~ 

30 

IF 

IN 

.GT. 233) 

15C.1A0 

AC 

IF 

( N 

.CT. A20 ) 

GC TO 60 

| 

IF 

TF 

.GT. Al AT 

GC TO 50 

i 

IF 

IN 

.GT. 39 A ) 

170,160 

!5C 

IF 

IN 

.GT. A15 ) 

19C, 180 

60 

IF 

TFT 

oG T * 7 8 1 1 

“gt: iu 70 


IF 

(N 

.GT. A21) 

2 1C ,2CC 

7 C 

IF 

IN 

.GT. 782) 

230,220 

C 0 • 

TCTB 

. SEC. 


80 

VPAe(2l = -3.AE 

-A - 1.5E-6*T 
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GG TG 260 

C 8. TC 2C. sec. 

90 VP AE (2 ) = -6.7E-3 - 8. £-5*1 
GO TG 260 

C 20. TG 7C. SEC. 

IOC VPA8 (2 ) = 2.604E-3 + I-6.684E-4 ♦ 6.16E-6*T)*T 
GC TC 260 

C 7C. TC 112. SEC. 

HO VPAB (2 ) = -2.17E-2 + 1.1E-4*T 
GC TC 260 

C 112. TC 112.5 SEC. 

120 VPAEil) = -1.3E-3 
VPAB12) = 8.2E-2 
VPAE13) = 5 .05E— 2 - 6.25E-4*T 
GC TC 260 

C 112.5 TC 116 SEC. 

130 VPAeiii = -l.E-3 

VPAB ( 2 ) = .1305 - l.E-3*T 
VPA6I3) = 5.05E-2 - 6.25E-4*T 
GC TC 260 

C 116. TC 116.5 SEC. 

140 VPAEUJ = -.24112 ♦ 2.07E-3*T 
VPAB 12 ) = 4.1441 — 3.56£-2*T 
VPAE13) = -5.29304 ♦ 4.544E-2*T 
GG TO 260 

C 116.5 TC 156.5 SEC. 

150 VPAEdJ = — 1.888125 E— 5 * 4.625E~7*T 
VPAB i 2 J * — 3. 8825E— 3 * 5.E-6*T 
VPAET3J = 8.365E-4 - l.E-6*T 
GG TC 260 

C 156.5 TC 157. SEC. 

160 VPABll) = 2.8368E-2 - 1.44E-4*T 
VPAB (2 ) = -1.1426 ♦ 5. 8E-3*T 
VPAB13) = .25216 - 1.28E-3*! 

GG TG 260 

C 157. TC 207. SEC. 

1 70 VP AB (2 ) * C. 

GG TO 240 

C 207. TC 207.5 SEC. 

180 VPAB <2 ) = -11.178 ♦ 5.4E-2 + T 
VPAB (3) = -.828 ♦ 4.E-3*T 
GG TC 260 
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C 2C7-5 TG 2 1C . SEC. 

190 VPABI2) = 2.7E-2 

VPAB 13 ) = -B.1E-2 ♦ 4.E-4*T 
GC TO 260 

C 210. TO 2T0. 5 SEC. 

200 VPA842) * 12.207 - 5.BE-2*T 

VPAB (3 ) = 1.011 - 4.8E-3*T 
GG TG 260 

C 210.5 TG 390.5 SEC. 

210 VPAB 1 2 ) = -2. £-3 

VPAB (37* 6.E-4 

GG TG 260 

C 390.5 TG 391. SEC. 

220 VPAB (2) = -1.9E-2 
VPAB (3) = 7.2E-3 
GG TC 260 

C 391. TG 400. SEC. 

230 VPAe(2) * 0. 

C ZERO CCNPGNENTS 
240 VPA6I3 ) - 0. 

250 VPABU) = 0. 

C BODY ACCELERATION 

260 if in .er 2 tn go to 290 — 

IF (N .GT. 140) GG TG 280 

IF IN .GT. 100) GO TG 270 

IF I N .GY. 80) 320.310 

270 IF (N .GT. 112) 340.330 

280 IF IN .GT. 232) 60 TC 370 

TF IN .GT.ZZ4TT60 ,350 
290 IF (N .GT. 414) GG TG 300 

IF IN .GT. 395) GO TO 400 

IF IN .GT. 394) 390,380 
300 IF IN .GT. 720) GG TC 430 
IF IN .GT. 415) 42C,41C 
CO. TO 40. sect 
310 APLBIl) * 14. ♦ .054* T 

APLBI3) = -.25 - .CC625*T 
GO TO 440 
C 40. TO 50. SEC. 

320 APieil) = 16.16 

APLBTJ) * -.5 

GO TO 440 
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C6*T 


C 50. TO 56. SEC. 

330 APLBtl) = 20.16 - . 

GO TC 44C 

C 56. TO 7C. SEC. 

340 APLBtl) = 10.08 ♦ .1*T 
GC TO 440 

C 70. TO 112. SEC. 

350 APLBtl) = 4.48 + .18*T 
GO TO 44C 

C 112. TC 116. SEC. 

360 APLBtl) » 2.24 ♦ .2*T 

A P L B 1 3 ) = 41.5 - • 3 75*T 
GO TC 440 

C 116. TO 136. SEC. 

370 APLBtl) * 2.24 + .2«T 

APLB ( 3 ) = -.956 - .CC9*T 
GC TC 440 

C 136. TO 197. SEC. 

380 APLBtl) = 29.44 

APLB t3 ) = .54 - .C2*T 
GC TO 440 

C 197. TO 197.5 SEC. 

390 APLBtl) « 11636.68 - 58.92*1 
APLE 1 3 ) = -1319.36 ♦ 6.68*T 
GC TO 440 

C 197.5 TO 2C7. SEC. 

400 APLBtl) = -.02 
APLB (3 ) = -.06 
GO TO 440 

C 2G7. TO 207.5 SEC. 

410 APLBtl) = -6106.52 + 29.5*T 
APL813) = -24.9 * .12*T 
GC TC 44C 

C 2C7.5 TC 36C. SEC. 

420 APLBtl) = i./t. 11402813064487 - 2.22358882 17603E-4*T ) 
APLBT3) = 0. 

GC TO 440 

C 360. TO 400. SEC. 

430 APLeti) = 29.43 

440 RETURN 
END 
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A. 3 OUTPUTS 


This subsection presents the outputs of a typical run. 


151 



FAILS I M _. . ... _ ..... 

02/24/73. C9.24.12. 

PROGRAM 'CONTROL 

OT 1 . 25000006-01 I ROM 2517566905 

N3MAX 2 NAM AX _ 4 . 

VEHICLE VIBRATORY MOTIONS 

AVBS 1 . OOOOOOOE +0 0 _1 . 0000 300S+00 1 . OOOOOOOE *00 

TAUV 4.00000000-02''" 4^00000000-02 4. OOOOOOOE -02 

KV 7.CC0COOGF+OO 7. OOOOOOOt+OO 7.0000000E+09 
KP 2. 5000000E + 01 2 . 50 0 JOCOE +0 1 2.5000000E+01 

RADIUS 6.C00C000E+00 


LMGNTt 0 LPL flT 0 N1>«AX I N2MAX 400 


9 . 8100000 E- 31 9 .jB 10000 OE-_pl 9.3 lODOODE-OJ. 

4 . 0000000 E-02 4. OOOOOOOE -02 4. OOOOOOOE -02 

7.0000000E-01 7 .0 00 0000 E — 01 7. OOOOOOOE -0 1 
2. 51000000 E-01 2 .5 00000 0_E-01_ 2_. 3 OOOOOOt-0 1 


Ui 

ts) 


INSTRUMENT 

ORIENTATION 

IRl 

1 

ThTRlI 

I R 2 

2 

THTK2I 

IP. 3 

3 

THTR21 

THTAI 

0 . 


THTGI 

0 . 




UNFAILED ERROR S 

SIGUFEA 8 • 6000000F + 0 l 1 .0 JOOOOOE+OO 7. OOOOOOOE *0 1 5. OOOOOOOE ♦01 5 .C OOOOOOfc *01 

SIGUFEG 8.3000000E-02 1 . 0003000F-02 1 . OOOOOOOE ■*■02 1 . OOOOOOOE +01 1 .OOOOOOOF «-01 2. OOOOOOOE-Ol 

2 .GO 3COOOF-02 2 .OOOOOOOE— 01 4. OOOOOOOE-02 

GACCI 2 • CCOCOOOE+OO GGYRGl 5 .OOOOOOOE+OO 


FAILURES 

NFAIL 

160C 

160C 

C 

C 

C 

c 

c 

0 

c 

c 

0 

c 


NINST 

1 

6 

0 

0 

0 

0 

0 

9 

0 

0 

0 

c 


MODt FAIL11 . FA1L2I 

1 _ 1 . OOOOOOOE +01 0. 

"l "" 1. OOOOOOOE ♦01 ~ 0. 

0 P. 0. 

0 G. 0. 

""o “ " "0. ” ■ " 0 . 

0 C. 0. 

0 0. 0. 

0 0 . o'. 

0 0. 0. 

0 0. 0. 

b r o. o. 

0 0. 0. 
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o. a - o .o. - a. 

C 0 0 0. 0. 

NOMINAL ALGORITHM. _ . . 

LALGO 1 

ADAPTIVE ALGORITHM 

LALG1 0 

FIETEENJiHRtSHCLU algorithm _ _ ] 

LALG2 f ' 


THAI 

I .99500005+03 

THG1 

3.01000005+00 

TAUPAQ 

1. OGOOOOOE+ 01 

TAUPGQ 

1.00000C-0E+01 

SQUARED 

LALG3 

LHEAL 

ERROR ALGORITHM 

I 

0 






| 

1 AK1I 

GK 1 1 
TAUPAQ 

2.2400000E+03 
4 .460 COG OE + 00 
l.CCOOOOOE+OI 

AK2I 
GK 2 1 
TAUPGQ 

4. 4444444E-01 
4 .44444445 -0 1 
1 .OOOOOOOE +01 

AK3I 
GK 3 1 
AKGOI 

1.9950000E+03 
7 . 9400000E+ GO 

C. 

AK4 I 
GK4I 

3.8095238E-01 

3.8095238E-01 

BAYESIAN 

LALG4 

DEC IS ICN THEORY 

l 

ALGORITHM 






j THAU 

TAUPAQ 

3 . 16000006+03 
l.COOGOOOE+Ol 

THA2I 

TAUPGQ 

3.I600000E+03 
I. 00000 00 £+01 

THUII 

6.3100000E+00 

THG2I 

1.2590000E+01 

MAXIMUM 

LALG5 

L IKEL 1 HOOD ALGORITHM 
0 






I MINIMAX 
! LALG6 
, LSQA 

ALGORITHM 

1 

1 LSQG 

1 






THAMI 

TAUPAQ 

9.780000QE+02 

l.OOOGOOOE+Ol 

ThASI 

TAUPGQ 

1 .OOOOOOOE +02 
I. 00000005+01 

THGMI 

1. 9530000E+00 

THGSI 

l.GOOOOODE-31 


| sequential algorithm 
LALG7 0 


I 

i 

1 
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| TIME 200.000 GYRO NO. 1 FAILED IN MODE 1 

! TIME 200.000 GYRO NO. 6 FAILED IN MODE 1 

! T IME 203.875 ALG. N O. 4 GYRO NO. 1 TKP E ALARM. 

priME 206.250 ALG. NO. 6 GYRO NO. L TRUE ALARM 


| TIME 220.375 ALG. NO. 6 GYRO NO. 6 TRUE ALARM 

j l TIME 343.875 ALG. NO . 3 GYRO NO. 6 TRUE ALARM 

j TIME 400.000 ALG. NO. 2 GYRO NO. 1 MISSED ALARM 

! TIME 400. 030 ALG. NO. 2 GYRO NO. 6 MISSED ALARM 


TIME 400. OOC ALG. NO. 3 GYRO NO. 1 MISSED ALARM 
TIME 400.000 ALG. NO. 4 GYRO NO. 6 MISSED ALARM 


i j TRIAL NO. 1. MAGNITUD ES O F ANGLE A ND VELOCITY ERRORS 

c n o o A n r\ n -j i. r\ n~ra a 1 i 


508.63 

0.00 

972.60 

978.51 

1039.77 

/ 

0.00 

943.32 

777.40 

1319.35 

13.0010 

7946.3733 

15.8938 

15.8011 

6.5079 


7951.6343 

11.9667 

10.7878 

12.5036 

! FDIC TIMES FOR 

TRIAL NO. 

1 



j 

0.000 

8.135 

5.704 

6.842 

1 

0.000 

.004 

. 006 

.007 


if IROM = 10975562961 


TOTAL TIME = 


132.147 





0.00 

425.02 

0 . GO 

G.00 

0.00 

119.96 

0.00 

0.00 

7946.3783 

11.7081 

7946.3783 

75»4Lii.7J>3_. 

7951.6343 

2.3058 

7951.6343 

7951.6343 


0.000 

13.692 

0.000 

0.000 

0.000 

.009 

0.000 

0.000 
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I 


MEAN. MAGNITUDES OF ANGLE 
*>08.83 0.00 

AN 

0 VELOCITY 

ERRORS 









. 


972.60 

973.51 

1939 

.77 


0 

.00 

“♦25.02 


s.oc- 

0.00 J 

0.00 


943.32 

777.40 

1319 

• 35 


0 

.00 

119.96 


0.00 

0.00 1 

13.0010 7946.3733 


13.3938 

15.3011 

8.5079 


7946.3 

7b3 

11.7081 

7 94 1. 

.37 33 

7946.3783 ( 

7951.6343 


11. 966 7 

10.787b 

12.5 >36 


7951.6343 

2.3058 

7 951 

• 63 43 

7951.6343 [ 

FDIC TIMES ECU ALL TRIALS 












I 

... .. .1 

0.000 


8.135 

5.7 04 

6. 

842 


0. 

000 

13.692 


0.000 

0.000 

0.000 


.004 

.006 

• 

007 


0. 

000 

.009 


0*000 

0.000 . 1 



0 Y P( 

: ALARM STATISTICS 




A CO EL fcROMRT ER ALA3M STATISTICS 

ALGORITHM 

l 

2 3 

4 5 

6 7 


6 

1 

2 

3 4 

5 

6 7 

3 ! 

ALL CLEAR 

o 

o'" 0 

' 0 0 

0 0 


0 

0 

C "" 

d 0 

0 

0 " "o 

" o ! 

CORRECTED FALSE ALARM 

0 

0 0 

0 0 

0 0 


0 

0 

0 

0 0 

0 

0 0 

o ! 

FALSE ALARM 

0 

0 0 

0 0 

_0 0 


0 

0 

,_0 

0 0 

0 

0 0 

- il i 

LOST ALARM 

0 

c 6 

0 0 

0 0 


o" 

0 

0 

0 0 

0 

0 0 

0 | 

TRUE ALARM 

0 

0 1 

1 0 

2 0 


0 

0 

0 

0 0 

0 

0 0 

' 0 i 

FORTUITOUS TRUE ALARM 

0 

0 0 

0 0 

U 0 


0 

0 

0 

0 0 

0 

0 0 

o ! 

FORTUITOUS ALL CLEAR 

0 

0 0 

0 0 

0 0 


0 

0 

0 

0 0 

0 

0 0 

i o 1 

MISSED ALARM 

0 

2 1 

l 0 

0 0 


0 

0 

0 

0 0 

0 

0 0 

’ 0 1 

MISSEO ALL CLEAR 

0 

C 0 

0 c 

0 0 

. 

0 

0 

0 

0 0 

0. _ 

0 0 

0 } 



APPENDIX B 

THE SEQUENTIAL ALGORITHM IN STEADY STATE 


The Kalman-Bucy filter equations for one of the three identical 
one-state filter equations derived in Paragraph 3. 1.2.4 are 

= x. + k.(y. - x i} (B-l) 


k. 

1 


m 


i 


m i + r 


p. = (1 - k.) xn i 



(B-2) 

(B-3) 

(B-4) 


m i+l = Pi + (B_5) 

where c|> , m, and q correspond to (j> 2 , m 2> and q 2 , respectively, in 
Eqs. 3-76, 3-77, and 3-78. 

If we combine Eq. B-l and B-4, we obtain 

J i+1 i i+ k i+1 (y i+1 -+ ij) (B-6) 

Taking the Z transform, and assuming k^ constant. 


* kz 

X ~ z - <(> (1 - k) ^ 


(B-7) 


The Z transform of a first-order constant coefficient filter of gain and 
time constant is 



(B-8) 
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Thus we see that, in the steady state 


k = K f [1 - e 


<t> (1 - k) = e 


JL 

T, 


_T_ 

T f 


(B-9) 


(B- 10) 


are the equations relating the steady- state Kalman filter to the constant 
coefficient filter. If we choose 4> as 


_T 

T 


4> = e 




(B - 1 1 ) 


then 




JL T 

T f" T * 


(B - 12) 


Now consider the covariance equation. From Eqs. B-2, B-3, and B-5, 


Pi+i x ( l - k i+i)( ,|,2 Pi + , >) 


(B-13) 


$ Pi + q 

e i+i "72 7 ~ 

<♦> Pj + q t r 


Using k i+ j from Eqs. B-14 and B-13 gives 

r(4> 2 Pi + q) 

P i+1 = ~2 7 7~ 

4> Pi + q + r 

From Eqs. B-14 and B-15, we have 


(B-14) 


(B-15) 


p i+l = r k 


i+1 


(B-16) 
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In steady state 


p. = p 


(B-17) 


and 


p = r k 


(B- 18) 


From Eqs. B-14, B-17, and B-18 we can obtain 


q 


[l - (1 - k) <t> 2 ] 
1 - k 


(B-19) 


Thus we can set up the Kalman-Bucy filter to have a desired steady- 
state time constant. First we choose the time constant and the 
measurement variance r. Then we select the state transition matrix by 
choosing (a typical value would be 2 t^). The steady- state Kalman-Bucy 
filter gain k is then given by Eq. B-12. The equivalent filter gain K f may 
then be found from Eq. B-9. We now find the initial state covariance p Q 
from Eq. B-18 and the state noise covariance q from Eq. B-19. The 
filter will now start and remain in steady state. 


If we wish to, we may now assign some other value to p Q (without 
changing q). The filter will then undergo a transient, eventually reaching 
the same steady state characterized by and 
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